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PREFACE 


This  report  is  the  final  report  on  the  tasking  by  Daniel  H.  Wagner  Associates  to  the 
Defense  Advanced  Research  Projects  Agency  (DARPA)  under  contract  MDA  972-9  l-C-0077. 
The  aim  of  the  research  is  to  develop  low  frequency  active/low  probability  of  intercept  acoustic- 
signals  and  detection  mechanism.  Our  approach  includes  implementing  the  discrete  wavelet  theory 
for  detection  and  classification  purposes. 

The  principal  investigator  is  Dr.  Robert  J.  Lipshutz.  Contributors  to  this  report  include  Dr. 
Eric  Dew  and  Mr.  Darren  Miner.  Mr.  Derek  Bernhart,  Ms.  Kirsten  Holmquist,  Mr.  Jason 
Kastner,  and  Ms.  Angela  Schillace  provided  additional  support  for  software  development,  graphics 
production,  computation  results,  and  the  write-up  of  the  User's  Guide  for  the  transmission 
simulation  system.  Consulting  services  were  provided  by  Dr.  Jay  Epperson  and  Dr.  Michael 
Frazier. 


SUMMARY 


This  final  report  presents  the  methodology  used  and  the  results  achieved  in  the  Low 
Frequency  Active/Low  Probability  of  Intercept  research  project  for  DARPA. 

THE  PROBLEM.  The  general  objective  of  this  work  is  to  develop  a  fast,  accurate  signal 
detection  methodology  based  on  Discrete  Wavelet  Theory.  The  questions  addressed  in  this  report 
include: 


•  Sonar  signal  choice.  We  discuss  the  general  types  of  sonar  signals  currently  in  use. 
The  choice  of  sonar  signals  must  take  into  account  Doppler  and  other  transmission 
effects,  probability  of  intercept,  and  type  of  information  one  wants  the  signal  to  relay. 
Physical  constraints  such  as  projector  characteristics  will  also  play  a  role  in  determining 
the  type  of  signals  one  can  send. 

•  Signal  detection.  One  needs  a  quantitative  measure  to  determine  whether  a  return 
signal  is  present  in  the  incoming  sample.  This  measurement  is  made  difficult  by  the 
presence  of  stochastic  noise  in  the  return  signal.  An  additional  complication  arises  in 
shallow  water  scenarios.  Here,  the  sound  channel  includes  multiple  paths  due  to  bottom 
bounce.  Multiple  paths  change  the  structure  of  the  signal  quite  significantly,  and 
sometimes  by  so  much  so  that  one  cannot  recognize  the  relation  between  the  return 
signal  and  the  prototype  signal. 

In  the  processing  end,  one  needs  to  know  where  the  signal  lies  within  the  detection 
window.  Such  information  yields  understanding  about  the  target's  position  or 
momentum. 

•  Choosing  a  detector.  Although  we  have  settled  on  the  optimal  Neyman-Pearson 
detection  algorithm,  its  implementation  in  the  compressed  version  is  basis  dependent. 
For  a  given  sonar  signal  choice,  one  would  want  a  basis  which  is  capable  of  efficiently 
detecting  the  return  of  this  signal. 

•  Algorithm  speed.  The  computation  must  be  fast  c. iough  to  handle  a  constant  stream 
of  incoming  information  in  real  time. 

METHODOLOGICAL  RESULTS.  Some  of  the  methodological  results  are: 
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•  Construction  of  a  multi-level  wavelet  basis  from  two  generating 
functions.  Examples  include  the  real  Shannon  basis  and  the  Daubechies  D4  basis. 
We  also  include  the  Haar  basis  for  mathematical  illustration. 

•  Construction  of  a  multi-level  wavelet  basis  in  which  the  p^-level 
wavelet  function  approximates  a  given  signal.  This  technique  allows  one  to 
tailor  the  wavelet  transform  to  detect  a  given  signal. 

•  Construction  of  a  multi-level  wavelet  basis  in  which  ail  the  wavelet 
functions  from  each  level  are  close  to  a  given  signal.  This  technique  is 
designed  to  yield  small  compression  loss  for  reproduction  and  detection  purposes. 

•  Construction  of  a  multi-level  wavelet  basis  in  which  wavelet  functions 
are  close  to  a  given  signal  and  to  all  rotations  of  it.  This  technique  allows 
one  to  tailor  the  wavelet  transform  to  detection  of  a  given  signal  in  a  robust  manner, 
allowing  for  errors  in  positioning  the  signal  in  the  sample  window. 

•  Construction  of  specialized  signal  types.  These  are  standard  sonar  signals 
having  varying  properties.  These  properties  include  being  invariant  or  very  sensitive  to 
Doppler  shifts. 

•  Description  of  the  compressed  Neyman-Pearson  detector.  We  give  formulas 
for  compressing  the  standard  Neyman-Pearson  detector  using  wavelets  so  as  to 
increase  the  speed  of  the  algorithm  without  significantly  degrading  performance. 

•  Description  of  the  multi-channel  detector.  This  detector  is  based  on  the 
compressed  Neyman-Pearson  detector  and  a  p-level  /2-closest  orthonormal  wavelet 
basis.  It  is  designed  to  detect  the  presence  of  any  translate  of  a  returned  signal  even 
when  there  is  Doppler  shifting. 

EMPIRICAL  RESULTS.  In  order  to  test  the  utility  of  wavelet-based  detectors,  we  created  a 
set  of  five  standard  signal  types  common  to  sonar  processing.  These  types  are  continuous  wave 
(CW),  hyperbolic  frequency  modulated  (HFM),  linear  frequency  modulated  (LFM),  pulse  (PUL), 
and  pseudo-random  noise  (PRN).  These  signals  are  depicted  in  Figure  V-2.  Individual  signals  of 
each  type  are  made  by  setting  certain  signal  parameters  like  the  bandwidth  or  the  sampling 
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frequency.  "Chip"  signals  are  constructed  by  concatenating  a  sequence  of  the  above  type  signals. 
Transmission  effects  are  simulated  on  a  chip  signal,  including  Doppler  shifting,  and  the  signal  is 
offset  to  model  the  effect  of  not  capturing  the  returned  signal  properly  in  the  receiving  sample 
window.  The  multi-channel  detector  methodology  uses  the  compressed  Neyman-Pearson  detector. 
The  basis  choice  used  in  this  detector  is  made  by  the  user.  We  show  that  orthonormal  /  -closest 
wave'et  bases  generally  outperform  the  more  standard  bases  like  the  Fourier  basis  or  the  delta 
(time)  basis.  Details  of  the  test  results  are  given  in  Chapter  VII.  The  major  empirical  results  are: 

•  The  /2-closest  wavelet  basis  outperforms  the  other  wavelet  bases  tested 
using  the  compressed  Neyman-Pearson  detector. 

•  The  /2-cIosest  wavelet  basis  performs  best  for  an  unaffected  signal  with 
no  offset  or  Doppler  shifting. 

•  Our  experiments  with  the  multi-channel  detector  were  successful  in  that 
the  detector  carried  out  all  the  tasks  requested.  However,  the 
experimental  results  were  not  enlightening.  The  rather  uninformative 
results  were  due  mainly  to  the  choice  of  the  prototype  signal.  On  the 
other  hand,  the  multi-channel  detector  showed  that  only  by  using 
wavelets  could  such  results  be  achieved.  Wavelets  allowed  for  compact 
storage  and  efficient  recursive  algorithms,  which  were  all  exploited  in 
building  the  multi-channel  detector. 

•  We  have  built  several  time-frequency  distributions,  or  spectrograms. 
These  spectrograms  give  us  a  new  way  to  analyze  signals.  The  wavelet 
spectrogram  is  another  way  to  see  the  power  of  the  time-frequency 
localization  of  wavelet  theory. 

FURTHER  RESEARCH.  Recommendations  for  further  research  include  the  following: 

•  Conduct  simulation  studies  using  real  world  noise. 

•  Conduct  simulation  studies  to  model  real  world  transmission  channels. 

•  Conduct  at-sea  tests  using  proposed  signals  and  analyzing  methodologies. 


Continue  study  on  other  potential  uses  of  these  /2-cIosest  wavelet  basis  for  signal 
processing. 

Continue  study  on  other  types  of  discrete  wavelets,  perhaps  those  that  are  more 
effective  for  time-frequency  analysis. 

Enhance  and  refine  the  current  transmission  models  used  in  the  simulation  software. 


CHAPTER  I 
INTRODUCTION 


This  document  is  the  final  report  on  the  Low  Frequency  Active  /  Low  Probability  of 
Intercept  Sonar  Signals  research  work  for  DARPA.  This  report  highlights  results  using  discrete 
wavelet  theory  in  the  development  of  near-optimal  detection  schemes.  We  also  detail  the 
construction  of  the  sonar  signals,  their  characteristics,  and  their  effectiveness.  Since  we  wish  to 
emphasize  practical  results  over  theoretical  results  in  this  report,  we  give  only  a  summary  of  the 
theory.  Our  main  premise  is  that  wavelet  theory,  especially  the  development  of  the  "/^-closest" 
discrete  wavelet  bases,  shows  substantial  performance  improvement  over  standard  bases  (i.e.,  the 
delta  basis  and  the  Fourier  basis)  for  signal  detection  and  computational  efficiency. 

Additionally,  we  have  created  a  simulation  system  to  simulate  the  effects  of  transmitting  a 
sonar  signal  underwater  -  including  Doppler  effects  and  bottom  bounce  -  to  test  the  effectiveness 
of  the  detection  method  as  well  as  to  highlight  the  characteristics  of  the  various  sonar  signals  that 
we  produce.  User  directions  and  details  of  the  front-end  are  given  in  the  User's  Guide:  Low 
Frequency  Active  Simulation  System,  Reference  [a]. 

1.  Overview  of  the  Contents 

A  complete  understanding  of  sonar  system  technology  is  beyond  the  scope  of  any  one 
report.  A  sonar  system  entails  everything  from  transmission  to  reception.  On  the  transmission 
end,  this  includes  the  hardware  choice  for  the  transmitter,  signal  choice,  power  capabilities, 
transmitter  characteristics,  and  signal  strength.  After  transmission,  one  has  to  deal  with  media 
effects  (i.e.,  effects  on  the  signal  during  its  trip  from  transmitter  to  receiver),  such  as  reverberation, 
Doppler  effects,  transmission  loss,  ambient  noise,  reflection  characteristics,  and  bottom  bounce. 
Finally,  at  the  receiving  end,  one  has  to  deal  with  signal  detection,  classification  of  signal  type, 
interpretation  of  the  return  signal  (e.g.,  has  the  signal  been  Doppler-shifted?),  receiver 
characteristics,  and  other  associated  problems.  We  should  add  that  we  are  considering  only  active 
sonar  signals,  that  is,  signals  that  the  receiving  team  has  sent  out.  Passive  sonar  technology  is  -  by 
default  -  simpler  at  the  transmitting  end,  but  much  more  difficult  at  the  receiving  end. 

Our  research  concentrates  on  two  main  aspects  of  the  sonar  system.  One  aspect  is  the 
signal  processing  phase  of  the  system  -  that  is,  the  receiving  end.  and  in  particular,  signal 
detection.  We  assume  that  the  ''received"  signal  -  one  which  is  ready  to  input  to  the  detector  -  has 


already  been  discretized  and  windowed.  Thus,  the  received  signal  is  represented  by  a  vector 

s  =  s( n |  =(s(0),s(l) . s (TV-  1)),  where  each  s(k)  is  a  complex  number.  Chapters  111  and  IV 

discuss  discrete  wavelet  theory  and  the  optimal  (and  near-optimal)  Neyman-Pearson  detection 
method,  respectively. 

The  second  aspect  deals  with  a  part  of  the  transmitting  phase  of  the  system,  namely  the 
signal  choice.  Here,  we  construct  and  evaluate  several  standard  signal  types,  generieally  called 
"chirps”  or  "chips."  Different  signal  types  have  different  characteristics,  some  of  which  are  useful 
for  certain  purposes  and  some  of  which  are  totally  useless  for  those  same  purposes.  For  example, 
some  signals  are  capable  of  detecting  Doppler  shifts:  some  are  invariant  under  Doppler  shifts. 
Thus,  the  latter  would  not  be  a  wise  choice  if  one  wants  to  know  whether  the  target  is  moving  or 
stationary.  These  topics  are  discussed  in  Chapter  V. 

Although  we  did  not  research  transmission  effects  (e.g.,  reverberation  and  bottom  bounce) 
in  the  same  vein  as  the  above  two  aspects,  we  did  develop  a  simulation  system  to  model  these 
effects.  Additionally,  we  developed  a  multi-channel  detector,  using  the  compressed  Neyman- 
Pearson  (near-optimal)  detector  technology,  to  simultaneously  detect  a  range  of  possible  Doppler 
shifts  of  an  active  signal.  Both  the  simulation  model  and  the  multi-channel  detector  are  detailed  in 
Chapter  VI.  Chapter  VII  gives  the  results  of  our  research  and  includes  many  graphs  supporting 
the  claim  that  certain  wavelet  bases  perform  better  than  other  bases  for  detection  purposes. 

The  areas  that  are  not  addressed  in  this  report  and  that  were  not  considered  during  the 
research  consist  mainly  of  hardware  and  software-hardware  interface  issues  (such  as  size  of 
window,  type  of  windowing  function,  and  sampling  rate).  However,  we  do  address  issues 
concerning  computing  time  and  memory  storage  capabilities.  We  strive  to  minimize  computing 
time,  and  this  is  achieved  to  some  extent  through  the  use  of  wavelets.  The  recursive  nature  of 
wavelet  functions  also  allows  one  to  minimize  storage  problems.  We  add  here  that  these  issues  are 
not  separate  from  the  above  issues  of  signal  detection  and  generation  and  transmission  simulation. 
Our  research  in  detection  methods  was  guided  in  part  by  selecting  those  algorithms  that  minimize 
computing  time.  The  multi-channel  detector,  on  the  other  hand,  was  created  only  after  realizing  the 
potential  storage  savings  that  wavelets  afford,  as  opposed  to  other  bases  like  the  Fourier  basis. 


2. 


Sonar  Signal  Technology 


Sonar  systems  are  prescribed  by  certain  parameters,  collectively  known  as  the  sonar 
parameters.  These  basic  parameters  are  needed  when  designing  the  system.  Since  we  are  not 
designing  the  system,  we  have  assumed  reasonable  values  for  these  parameters.  Some  parameters 
are  unnecessary  for  us  to  specify;  others  are  specified,  but  need  not  be  strictly  adhered  to.  The 
classes  of  sonar  parameters  considered  are: 

•  Projector  size.  This  measures  the  power  output  capabilities  of  the  sonar  projector.  This 
parameter  is  used  mainly  in  the  simulation  model. 

•  Projector  characteristics.  This  describes  the  frequency  spectrum  of  the  projector.  This 
parameter  is  also  used  in  the  simulation  model. 

•  Transmission  effects.  These  parameters  include: 

-  Bottom  bounce  and  reverberation  (for  shallow  water  scenarios). 

-  Doppler  effects.  We  assume  relative  speeds  <  15  knots  (kts). 

-  Distance  attenuation.  Total  source-to-target-to-receiver  distances  are  in  the  50-100 
nautical  mile  (nm)  range. 

-  Noise.  There  are  a  number  of  noise  sources.  One  is  biologies,  like  schools  of  fish 
or  whales.  Shipping  also  produces  noise,  mainly  from  the  blades  and  the  trailing 
eddies.  Natural  convection  of  fluids  also  produces  a  characteristic  noise. 
Reverberation  of  all  the  above  noise  upon  topography  (e.g..  rocks)  and  biologies  (a 
school  of  fish  can  reflect  noise)  also  produces  noise. 

The  noises  described  above  are  not  entirely  stochastic:  that  is.  each  type  is 
characteristic  of  the  source.  For  example,  propellor  noise  is  distinctly  in  the  50  or 
60  Hertz  (Hz)  range.  There  is  another  source  of  noise  that  is  purely  stochastic. 
This  source  is  the  ambient  "white"  noise  that  is  present  in  every  system  (including 
the  machines  that  transmit,  receive,  and  analyze  signals).  In  the  detection 
algorithm,  we  assume  only  white  noise.  We  have  not  modeled  real  noise  for  use  in 
the  detection  algorithm,  although  real  noise  can  be  applied  in  the  simulation  system. 
Such  implementation  requires  further  research  into  stochastic  noise  modeling. 


•  Target  strength  and  reflective  properties.  This  measures  the  effective  signal  reflection  of 
the  target  based  on  its  orientation  relative  to  the  signal  path  and  on  the  reflective  surface 
characteristics  of  the  target. 

3.  Broad  Technical  Issues 

As  stated  in  Section  1  above,  we  concentrate  on  two  main  aspects  of  sonar  systems 
technology.  The  first  is  detection,  and  the  second  is  signal  creation.  The  primary  issue  concerning 
detection  is  the  following  question: 

Is  there  a  signal?  Even  in  an  active  sonar  system,  in  which  the  time  of  transmission  is 
known,  it  is  impossible  to  know  exactly  when  the  signal  will  return.  Thus,  the  receiver  is 
constantly  receiving  something  that  we  generically  call  "incoming  signal."  Embedded  in  this 
incoming  signal  may  be  the  active  signal  that  was  transmitted.  One  only  has  to  detect  it  above  the 
ever-present  noise.  However,  the  active  signal  has  undergone  significant  changes  due  to  the 
transmission  effects.  (By  modeling  the  transmission  effects,  we  hope  to  give  an  accurate 
simulation  of  a  returned  signal  to  test  the  detector.) 

The  detector  is  a  variation  of  the  Neyman-Pearson  detector.  The  Neyman-Pearson  detector 
is  optimal  in  the  sense  that  for  a  fixed  probability  of  false  alarm  and  Gaussian  noise,  the  detector 
has  the  highest  probability  of  detection.  The  near-optimal  versions,  which  we  use,  sacrifice 
optimality  for  computing  speed.  Nevertheless,  our  research  shows  that  wavelet-based  compressed 
Neyman-Pearson  detectors  do  not  lose  much  optimality  in  return  for  substantial  gains  in 
computing  speed. 

A  subtler  issue  is  the  location  of  the  signal  within  the  window.  The  location  of  the  active 
signal  within  the  window  yields  range  information,  which  is  essential.  On  the  other  hand,  the 
detector  should  be  location-invariant.  This  means  the  detector  should  detect  the  presence  of  the 
active  signal  positioned  anywhere  in  the  window.  Certain  bases  (such  as  the  delta  basis)  give 
better  range  information  than  other  bases,  but  the  drawback  is  that  detectors  using  this  delta  basis 
do  not  perform  equally  well  when  the  signal  is  translated  within  the  window. 

This  trade-off  is  part  of  a  broader  constraint  known  as  the  Uncertainty  Principle.  The 
Uncertainty  Principle  forbids  the  simultaneous  determination  of  both  range  and  spec.:  information 
with  complete  precision.  To  deal  with  this  constraint,  one  devises  active  signals  that  can  yield  the 
desired  information  on  range  or  speed.  In  our  situation,  the  detector  is  built  upon  the  basis  chosen 
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to  represent  the  vector  space  (this  is  explained  in  Chapter  IV).  and  different  bases  yield  different 
information.  Wavelet  bases,  in  general,  are  characterized  by  their  ability  to  yield  both  range  and 
directional  information,  albeit  the  Uncertainty  Principle  forces  the  information  to  be  imprecise. 

We  create  several  types  of  sonar  signals  and  demonstrate  their  characteristics,  especially 
whether  they  yield  range  or  direction  information.  These  characteristics  are  shown  by  a  surface 
plot  called  a  Doppler  Ambiguity  display.  Such  plots  (see  the  figures  at  the  end  of  Chapter  V) 
show  how  "different"  a  signal  is  when  it  is  shifted  in  time  (thus  yielding  range  information)  or  in 
frequency  (thus  yielding  Doppler  shift,  or  direction,  information). 

One  of  the  early  directions  of  this  research  was  to  use  wavelets  to  generate  active  signals 
with  desirable  characteristics.  Subsequent  research  indicated  that  wavelets  would  provide  no 
significant  advance  in  signal  design.  Nonetheless,  we  devised  a  method  to  use  any  signal  to 
generate  wavelet  bases  for  detection  purposes.  This  is  described  in  Chapter  III. 

4.  Chapter  Summary 

We  start  with  a  brief  mathematical  and  historical  background  in  the  following  chapter. 
Chapter  III  gives  our  main  mathematical  results  in  discrete  wavelet  theory.  This  chapter,  as  well  as 
this  whole  report,  will  emphasize  the  results  and  will  delve  into  the  mathematics  only  to  give  the 
motivation  for  the  procedures.  Chapter  IV  discusses  the  detection  algorithm  and  includes 
quantitative  results  of  the  efficacy  of  certain  discrete  wavelet  bases.  Standard  sonar  signals,  along 
with  their  characteristics,  are  described  in  Chapter  V.  Chapter  VI  discusses  the  transmission 
simulation  model.  Chapter  VII  includes  an  explicit  example  of  its  use  in  conjunction  with  various 
detector  bases  as  well  as  with  the  multi-channel  detection  method.  Chapter  VII  consists  of  the 
results  of  the  experiments  and  tests  described  in  Chapter  VI.  We  end  with  a  discussion  of  possible 
future  research  directions  in  Chapter  VIII. 
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CHAPTER  II 


BACKGROUND  AND  NOTATIONS 


This  chapter  provides  the  motivation  and  notation  used  throughout  the  rest  of  this  report.  In 
Section  1,  we  trace  the  evolution  of  a  signal,  from  reception  through  the  analyzer  right  up  to  the 
detector.  Section  2  discusses  -  in  more  detail  -  the  mathematics  of  discrete  signal  processing.  The 
details  are  essential  for  the  notations  and  definitions  used  in  other  sections  of  the  report.  Section  3 
presents  a  brief  overview  of  the  historical  results  in  classical  harmonic  analysis,  which  led  to  the 
development  of  modern-day  signal  processing.  We  end  this  chapter,  in  Section  4,  with  a  brief 
historical  introduction  to  wavelets  and  its  relation  to  signal  processing. 

1.  Signal  Analysis 


A  signal  is  a  real-valued  function  of  one  variable,  usually  the  time  variable.  For  most 
purposes  we  do  not  have  to  assume  this  function  is  smooth  or  even  continuous,  but  we  will 
assume  the  signal  has  finite  energy.  Thus,  we  consider  the  signal  as  a  function  s(t)  e  L2(R)  (the 
function  space,  L2(R),  is  a  Hilbert  space)  that  is,  we  require 


|s(r)| 2  d  t<o°. 


As  briefly  discussed  in  Chapter  I,  this  signal  is  our  "incoming"  signal;  it  is  what  the  receiver  picks 
up  at  all  times.  In  this  case,  the  finite  energy  requirement  is  only  an  approximation,  since  at  all 
times  there  is  some  nontrivial  ambient  noise.  This  noise  is  frequently  above  some  threshold  so 
that  the  integral  over  all  time  would  indeed  be  infinite. 

Modern  signal  analysis  is  almost  exclusively  conducted  with  digital  computers.  This 
means  the  first  step  for  the  signal  after  reception  will  be  discretization.  For  this  process,  a 
"sampler"  picks  off  the  amplitude  of  the  signal  at  discrete  intervals.  We  are  not  concerned  with  the 
hardware  aspects  of  this  process  the  end  result  is  a  discrete-time  signal,  represented  by  an  infinite 
sequence,  s  =  s[n]  =(...,  s(— 2),  s(-l),  s(0),  s(l),  s(2), ...  )  (We  denote  the  vector,  s,  in  the 

standard  basis  by  sfn]  and  s(n)  to  denote  the  n*h  coordinate  of  s[n].) 

One  can  just  as  well  consider  this  infinite  sequence  as  a  real-valued  function  of  one  discrete 
variable.  In  the  continuous  case,  the  variable  is  a  real  variable.  In  the  discrete  case,  the  variable  is 
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an  integer  variable.  Thus,  the  signal  is  now  considered  as  a  function,  s[n]  €  /2(Z).  The  relation 
between  the  continuous  and  the  discrete  version  of  the  signal  is  given  by 


s[n]  —  s[nTsamp],  n  £  Z. 

The  sampling  period  is  Tsanip.  Its  inverse,  vsamp  =  VTsamp' is  the  sampling  frequency.  One  often 

fixes  the  sampling  frequency,  as  required  by  the  Nyquist  Sampling  Theorem,  to  be  greater  than 
twice  the  highest  expected  frequency,  which  can  be  determined  by  the  frequency  spectrum  of  the 
signal.  In  signal  processing  parlance,  s[n]  is  called  the  discrete-rime  signal. 

The  discrete-time  signal  is  now  ready  to  be  analyzed.  Analysis  involves  sending  the 
discrete-time  signal  through  a  filter1.  The  type  of  filters  we  employ  are  finite  impulse  response 
filters.  Finite  impulse  response  filters,  or  FIR  filters,  use  only  a  finite  sequence.  An  example  of  a 
FIR  filter  is  the  Discrete  Fourier  Transform  (DFT).  We  require  finite  sequences  because  we  do 
use  DFTs  in  our  computations  and  because  the  detector  cannot  process  an  infinite  stream  of  data  in 
real  time;  thus,  storage  problems  arise  as  the  detector  bottlenecks  the  whole  process. 

The  process  of  chopping  an  infinite  sequence  into  a  finite  sequence  is  called  windowing. 
This  entails  taking  a  windowing  function,  w[n]  e  /2{Z),  to  "mask"  a  section  of  the  discrete-time 
signal.  A  windowing  function,  or  window,  is  any  function  that  is  nonzero  in  a  finite  interval.  That 
is. 


,  f  *0forn=0,  1 . N-  1 

W  n  I  =  0  for  ail  other  n  e  Z  . 

We  assume  w(n)  >  0  for  n  =  0,  1 . N  -  1.  The  integer,  N,  is  called  the  window  size.  In  this 

report,  N  will  always  denote  the  window  size.  Also,  the  window  size  is  often  set  to  be  a  power  of 
2,  that  is,  N  =  2L,  for  algorithmic  speed  considerations.  Again,  N  is  usually  chosen  larger  than  the 

N_ 

sampling  frequency,  because  the  discrete  frequency  for  this  window  size  is  bounded  above  by 
thus  minimizing  information  loss.  In  the  research  conducted  for  this  report,  we  chose  the 
rectangular  windowing  function. 


,  ,  I  1  for  n  =  0,  l,...,yV-  1 
W(n)=i«  otherwise 


1  A  filter  is  a  generic  term  for  an  analyzing  processor  -  a  black  box.  so  to  speak.  Its  root  comes  from 

the  early  processors  which  literally  did  filter  out  the  high  or  the  low  frequencies. 
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We  chose  this  primarily  for  convenience.  To  make  the  discrete-time  signal  finite,  one  computes 
the  product  to  get  the  windowed  signal,  sw[n|,  where  sw(n)  =  w(n)  s(n)  for  all  n.  If  the  window  is 
the  rectangular,  then  there  is  no  need  to  compute  N  multiplications.  All  we  have  to  do  is  set 

sw(n)  =  s(n)  for  n  =  0 . N  -  1  and  set  sw(n)  =  0  otherwise.  Since  sw(n)  =0  for  n  *  0 . N  -  1 , 

we  consider  this  windowed  signal  as  a  finite  sequence,  which  we  also  label  as  s[n],  abusing  the 
notation  somewhat.  We  will  only  consider  finite  sequences  for  our  signals.  Such  signals  are 
called  discrete  signals.  A  source  for  discrete  signal  processing  is  Reference  [b]. 

2,  Some  Mathematics  of  Discrete  Signal  Processing 

The  discrete  signal,  s[n],  can  be  considered  in  a  number  of  ways.  One  way  is  to  consider 
s[nl  as  a  real-valued  function  of  one  integer  (modulo  AO  variable.  That  is,  s[n]  e  /2( Z/v)  where  Z/v 
denotes  the  integers  modulo  N.  Another  way  to  interpret  the  finite  sequence  s[n]  is  as  an  N- 
dimensional  vector.  In  this  case,  we  also  generalize  a  bit  by  viewing  s[n]  as  a  complex-valued 
vector.  This  generalization  facilitates  using  the  Fourier  transform,  which  is  complex-valued. 
Indeed,  we  could  have  started  at  the  very  beginning  and  considered  the  signal  s(t)  as  a  complex¬ 
valued  function  from  the  start.  We  choose  to  view  s[n]  as  a  vector  because  we  can  use  linear 
algebra  techniques  to  understand  the  problems.  However,  there  are  nonlinear  facets  involved  with 
discrete  signal  processing.  In  those  instances,  we  view  s[n]  as  a  complex-valued  function  on  the 
integers  0,  1, ...,  N  -  1. 

Now  that  we  have  established  s[n]  as  a  vector  in  complex  /V-dimensional  vector  space, 
denoted  by  V-  we  have  to  determine  the  signal  processing  equivalents  to  all  the  standard 
linear  algebra  objects,  like  bases,  matrices,  linear  operators,  and  so  on. 

Finite-dimensional  vector  spaces  are  best  viewed  independent  of  a  basis.  In  this  regard,  a 
signal  (i.e.,  vector)  s  no  longer  has  the  implication  of  a  function  of  time.  In  our  notation,  s  e  %  is 
the  signal  in  a  basis-free  notation.  The  standard  basis  for  the  vector  space  is  deemed  the  delta,  or 
time,  basis.  We  write  s[n]  e  %  to  represent  the  signal  in  the  standard  basis.  Since  all  signals  in 
time  are  real,  the  coordinates  of  s[n]  are  all  real-valued.  That  is, 
s[n]  =  (s(0),  s(  1 ),  s(2) s (N  -  1 )),  where  s(k)  €  R  for  all  k  =  0,  1, 2 N  -  1. 


X 


Another  important  basis  for  the  vector  space  is  the  Fourier  basis.  The  signal  in  the  Fourier 
basis  is  denoted  by  s[n].  The  relation  between  the  Fourier  basis  and  the  standard  basis  is  given  by 
the  Fourier  and  inverse  Fourier  transforms: 


N-  1 

s(n)  =  ^  s(m)  e-2  71 ' mn,N 

m  =  0 


and 


N  -  1 

s(n)  =  —  ^  s(m)  e2  n '  mn/A/. 

^m=0 


An  important  non-linear  algebra  operation  on  M)  is  the  convolution  of  two  signals.  For 
s[nj  and  t[n]  in  define  the  convolution,  (s*t)[n],  by 

N  -  1 

(s*t)(n)  =  s(k)  t(n  -  k), 

k  =  0 

where  the  index  is  computed  modulo  N.  The  main  convolution  result  is  the  convolution  theorem: 
Theorem  i:  (u*v)(n)  =  u(n)  v(n),  n  =  0,  1, ...  ,N~  1. 

See  Reference  [b|  for  a  proof.  The  multiplication  of  the  Fourier  transforms,  u|n]  and  v[n],  is 
coordinate-wise,  which  we  call  the  outer  product  of  the  two  vectors.  We  repeat  that  neither 
convolution  nor  the  outer  product  is  a  linear  operation. 

Additional  linear  algebra  facts  are  given  in  Reference  [c] .  We  now  summarize  the 
important  linear  operators  used  in  the  following  chapters.  Linear  operators  are  the  linear  algebra 
analogue  of  linear  systems  in  signal  processing.  Special  types  of  linear  systems,  such  as  linear 
time-invariant  systems,  are  just  special  types  of  linear  operators.  Some  systems  in  signal 
processing  are  not  linear.  These  correspond  to  nonlinear  operators,  and  in  this  case,  we  no  longer 
view  sfn]  as  a  vector;  we  view  it  as  a  function  whose  domain  consists  of  N  numbers. 
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Our  interpretation  of  linear  operators  here  differs  from  standard  linear  algebra  definitions  in 
one  very  crucial  aspect.  This  main  difference  is  that  our  operators  are  given  in  matrix  form  (no 
difference  here  with  most  linear  algebra  interpretations),  but  independent  of  basis.  In  this  respect, 
our  linear  operators  are  not  technically  linear  operators  in  the  linear  algebra  sense;  that  is,  we  name 
the  operator  by  what  it  does  to  a  vector  (expressed  in  any  basis)  not  by  what  it  does  to  a  vector 
(expressed  with  respect  to  a  given  basis).  This  is  a  very  important  difference.  The  reason  for  the 
difference  is  that  their  algorithmic  implementation  requires  us  to  do  so. 

The  first  linear  operator  is  the  rotation  by  k  operator.  This  operator  is  denoted  by  Rk.  It 
takes  a  vector  (expressed  in  whatever  basis)  and  translates  the  indices  by  k.  Specifically, 

(Rk  s)(n)  =  s(n  -  k). 


We  assume  here  that  the  basis  is  the  standard  basis,  which  is  where  this  operator  is  most  often 
used.  Signal  processing  engineers  usually  call  this  operator  the  delay  (by  k)  operator.  We  also  call 
this  operator  the  translation  (by  k)  operator.  Note  that  translation  by  N  is  the  identity  operator. 


Two  other  operators  coming  from  signal  processing  are  the  upsampling  and 
downsampling  operators.  They  are  denoted  by  U  and  D,  respectively,  and  are  defined  as  follows: 


(U  s)(n)  = 


s(y)  if  n  is  even 
0  if  n  is  odd 


,  n  =  0,  1,  2, ...,  2/V  -  1, 


and 


(D  s)(n)  =  s(2 n),  n  =  0,  1,2 . &-1. 

2 

Note  that  the  length  of  the  vector  doubles  when  upsampled  and  halves  when  downsampled. 
Again,  the  vector  is  expressed  here  in  the  delta  basis,  but  the  upsampling  and  downsampling 
operators  work  the  same  way  with  respect  to  any  basis.  Because  these  operators  take  vectors  from 
one  vector  space  to  a  different  vector  space,  one  must  be  very  careful  when  defining  the  bases  in 
the  respective  spaces. 

The  Fourier  transform  of  the  above  two  operators  results  in  the  spreading  and  folding 
operators.  These  are  denoted  by  S  and  F,  respectively,  and  are  defined  as  follows: 
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and 


for  n  e  {(),  1, N  -  1 } 

for  n  €  (A,  /V  +  1 . 2 N-  1}  , 


(S  v)(n)  = 


v(n) 

v(n  -  N ) 


(F  v)(n)  =  v(n)  +  v(n  +  N  ),  n  <=  { 0,  1 , ...,  N  ~  1 ) . 

Reference  [c]  details  the  relationships  between  the  upsampling  and  the  spreading  operators  and 
between  the  downsampling  and  the  folding  operators.  We  will  assume  knowledge  of  these 
relationships  in  the  following  chapters. 

3.  Harmonic  Analysis 

The  decomposition  of  functions  into  sines  and  cosines  as  described  by  J.  Fourier  at  the  end 
of  the  19  th  century  introduced  the  field  of  harmonic  analysis,  abstractly  culminating  with  the 
current  study  of  Hilbert  spaces.  On  the  practical  side,  there  are  a  number  of  reasons  for 
decomposing  a  function  -  or  in  our  case,  a  signal  -  into  trigonometric  components.  A  few  of 
them  are: 

1.  Frequency  Analysis.  One  would  like  to  know  the  tonals  making  up  the  signal. 

2.  Noise  discrimination.  With  a  priori  knowledge  of  the  signal's  frequency  range,  one 
may  zero  out  all  contributions  from  outside  the  range,  since  they  are  essentially  noise. 

3.  Data  storage  and  retrieval.  A  pure  tonal  has  infinite  support  in  the  time  domain. 
However,  its  Fourier  transform  has  compact  support.  Thus,  storage  is  simpler  in  the 
Fourier  domain. 

The  Fourier  transform  is  given  by  the  formula  (the  smooth  analogue  to  the  above  discrete 
Fourier  transform), 


f(co)  = 


f(r)  e-  m  dr. 


1 1 


although  the  integral  is  sometimes  normalized  by  scaling  with  The  inverse  transform  is  given 
by 


f (r)  =  —  I  f(to)  eic0'  do), 

2nJ 

J  .00 

where  the  scaling  constant  applies  when  the  above  Fourier  transform  is  used.  This  is  the 
continuous  analogue  to  the  discrete  Fourier  transform  given  in  Section  2. 

Current  applications-related  work  in  harmonic  analysis  include  time-frequency  (T-F) 
distributions.  T-F  Distributions  are  surface  plots  of  a  signal,  given  in  two  variables  corresponding 
to  time  and  frequency.  The  time-dependent  Fourier  transform,  or  the  windowed  Fourier 
spectrogram,  is  one  example  of  a  time-frequency  distribution.  Another  well-known  time- 
frequency  distribution  is  the  Wigner-Ville  distribution.  We  show  some  examples  of  these  T-F 
distributions  in  Chapter  V. 

In  our  research,  we  have  also  developed  a  wavelet-based  T-F  distribution.  This  distribution 
does  not  perform  as  well  as  the  windowed  Fourier  spectrogram,  but  it  is  much  faster,  since  it  does 
not  require  multiple  uses  of  the  Fourier  transform  over  a  number  of  windows.  Reference  fd]  gives 
a  thorough  treatment  of  the  various  T-F  distributions. 

Our  application  of  harmonic  analysis  comes  in  Chapter  V,  which  deals  with  signal 
generation.  Another  analyzing  device  we  use  is  the  Doppler  Ambiguity  Display.  This  is  also  a 
surface  plot,  much  like  the  T-F  distributions.  The  Doppler  Ambiguity  Display  is  not,  however,  a 
distribution  plot.  It  measures  how  "different"  a  signal  is  after  it  is  shifted  in  time  or  frequency  (by 
Doppler),  compared  to  its  unshifted  self.  The  difference  is  measured  by  computing  an  integral  of 
an  inner  product.  A  more  complete  explanation  is  given  in  Chapter  V,  as  well  as  in  Reference  [c]. 
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4. 


Wavelet  Theory 


One  reason  for  the  creation  of  T-F  distributions  is  to  study  nonstationary  signals. 
Stationary  signals  are  those  with  a  fixed  frequency  spectrum  over  all  time.  It  is  clear  that  most 
signals  (acoustic  signals,  in  particular)  are  nonstationary.  In  fact,  many  signals  are  transient. 
Finding  the  frequency  spectrum  of  a  nonstationary  signal,  then,  no  longer  has  its  original  meaning. 
Thus,  one  looks  at  the  windowed  Fourier  spectrogram,  which  assumes  the  signal  is  stationary  (i.e., 
has  a  fixed  spectrum)  over  a  small  period  of  time.  This  process  attempts  to  localize  (in  time)  the 
frequency  distribution,  because  the  Fourier  transform  (i.e.,  the  Fourier  basis)  has  no  time 
localization. 

What  is  needed,  then,  is  a  basis  which  has  time  localization.  Of  course,  the  delta,  or  time, 
basis  provides  perfect  time  localization.  However,  the  delta  basis  has  no  frequency  localization. 
This  means  the  delta  basis  offers  no  information  on  the  frequency  spectrum  of  the  signal,  which  is 
essential  if  one  were  to  try  to  quantitatively  distinguish  between  similar  looking,  but  different, 
signals.  For  example,  compare  a  linear  frequency-modulated  signal  to  a  hyperbolic  frequency- 
modulated  signal;  the  two  look  similar  in  time,  but  very  different  in  frequency  (see  Chapter  V  for 
illustrations). 

The  reason  for  these  constraints  is,  again,  the  Uncertainty  Principle.  The  Uncertainty 
Principle  is  quantitatively  given  by  the  inequality, 

afar  ^  nL 

where,  assuming  a  normalized  signal,  a(:  is  the  variance  in  the  frequency  v  and  and  of  is  the 
variance  in  the  time  t.  Since  al  =  0  in  the  Fourier  transform  of  a  signal  and  aj  =  0  in  the  time 
series  representation  of  a  signal,  the  complementary  variances  must  be  infinite. 

Thus,  the  inadequacies  of  both  the  time  domain  and  the  Fourier  (frequency)  domain  for 
signal  processing  led  researchers  to  consider  domains  (i.e.,  alternative  bases  for  L2(R))  that 
capture  a  bit  of  both.  The  Uncertainty  Principle,  people  recognize,  limits  any  domain  from 
precisely  determining  both  time  and  frequency  information.  In  this  regard,  there  is  no  optimal 
domain  in  which  to  perform  signal  analysis. 
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Results  in  this  area  led  to  the  study  of  wavelet  theory.  Initial  works  by  Gabor  based  on 
works  of  Wigner  (Reference  [e]),  and  many  others  considered  a  set  of  functions  in  time 
consisting  of  dilations  and  translations  of  a  single  scaling  function.  The  translates  would  operate  in 
the  same  vein  as  translates  of  a  delta  basis  function  (i.e.,  the  delta  function):  give  temporal 
information.  The  dilations  correspond  to  translations  in  the  Fourier  domain,  thus  yielding 
frequency  information.  This  set  of  functions  -  provided  the  functions  together  satisfy  some 
orthonormality  and  admissibility  conditions  -  forms  a  wavelet  basis  for  the  Hilbert  space. 

Many  researchers  in  wavelet  theory  currently  follow  the  Multiresolution  Analysis 
construction  developed  by  Mallat  and  furthered  by  Daubechies,  Mallat,  Meyer,  and  many  others, 
see  References  [f],  [g],  and  [h].  Here,  the  Hilbert  space  L2( R)  is  decomposed  into  a  nested 
sequence  of  closed  subspaces: 


•  •  •  V.i  c  V.\  c  T0  c  V\  c  V2  c  •  •  •  c  L2( R) 


such  that: 


(0)  =nvk=  lim  Vfc 

k— >-00 


and 

L2(R)  =  kjVk=  lim  Tk, 

k — >00 

where  the  overline  denotes  closure  (in  the  sense  of  Cauchy  sequences).  From  this,  there  exists  a 
unique  scaling  function  <])o(x)  e  To  such  that  for  <pk(x)  =  ^ <po(2kxl  ke  Z,  the  set 
^<t>k(x  -  2“  km)  I  m  e  Zj  spans  for  all  k. 

Define  complementary  subspaces  Wk  =  Tk  c  Tk+i>  k  e  Z,  so  that  VTk  ©  Tk  =  14+1.  Then, 
the  VV'k's  are  mutually  orthogonal,  and  the  direct  sum  of  them  all  spans  L2{ R).  Each  orthogonal 
subspace  is  in  turn  spanned  by  a  complementary  set  of  functions,  jvj/k(x  -  2~km)  |  m  e  zj. 

These  functions  are  called  wavelets.  The  V|/k's  are  constructed  to  complement  the  (j)k's.  See 
References  [f],  (g),  and  [i]  for  detailed  expositions  of  the  Multiresolution  Analysis. 
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CHAPTER  III 


DISCRETE  WAVELET  THEORY 

In  this  chapter,  we  cover  the  main  theoretical  results  of  discrete  wavelet  theory.  We  begin 
in  Section  1  with  a  brief  overview  of  discrete  wavelet  bases  and  how  they  are  constructed.  This 
description  includes  the  filter  bank  design  to  implement  the  wavelet  transform.  In  Section  2,  we 

-y 

study  the  construction  of  the  /  -closest  wavelet  and  its  corresponding  wavelet  basis.  Examples  in 
Chapter  V  show  the  utility  of  this  wavelet  basis  for  detection.  We  end  with  Section  3,  giving 
another  application  of  the  /  -closest  wavelet  basis:  computing  the  \jr-binary  decomposition  of  a 
signal.  This  decomposition  seems  to  be  useful  in  detecting  translates  of  a  given  signal.  There  are 
examples  of  its  use  in  Chapters  V  and  VII. 

1.1.  Discrete  Wavelet  Basis 

The  term  "wavelets"  generically  refers  to  a  set  of  functions  satisfying  certain  properties. 
These  properties  are  detailed,  for  the  continuous  case,  in  References  [f],  [g],  and  many  others. 
Some  of  these  properties  for  the  continuous  case  do  not  carry  over  to  the  discrete  case.  We  are 
concerned  specifically  with  the  discrete  case.  Nonetheless,  the  basic  premises  do  carry  over, 
namely  that  the  set  of  functions  forms  an  orthonormal  basis  and  that  translates  and  dilates  of  the 
functions  look  "alike." 

The  principal  difference  between  our  discrete  wavelet  functions  and  the  discrete  wavelet 
functions  from  some  of  the  literature  (see  References  [f],  [j],  [g])  is  that  our  discrete  wavelet 
functions  are  actually  discrete  functions  where  as  others'  discrete  wavelet  functions  are  just 
continuous  functions  that  have  been  discretely  sampled.  The  difference  is  that  we  can  guarantee 
normality  while  the  others'  wavelets  have  to  be  renormalized.  Furthermore,  our  version  allows  for 
direct  application  of  linear  algebra  and  is  easier  to  implement,  compared  to  the  other  types.  The 
drawback  is  the  lack  of  resolving  power.  The  discretely  sampled  continuous  wavelet  functions  can 
resolve  to  any  magnification.  Our  discrete  wavelet  basis  can  resolve  to  a  resolution  power  of,  at 
most,  log2  N  (which  is  10  for  N  =  1024). 

Some  of  the  wavelet  functions  are  illustrated  at  the  end  of  this  subsection.  Of  note  is  the 
Daubechies  D4  wavelet  basis.  We  took  the  algorithm  given  in  Reference  [g]  and  created  a 
"compactly  supported,"  orthonormal  wavelet  basis.  The  wavelet  basis  functions  are  compactly 
supported  in  the  continuous  sense:  as  functions  in  one  real  variable.  In  our  case,  all  our  wavelet 
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basis  functions  are  always  compactly  supported,  since  they  are  functions  in  one  integer  (modulo 
AO  variable.  From  the  pictures.  Figures  Ill- 1 ,  1II-2,  and  111-3,  we  see  that  both  the  real  Shannon 
wavelet  basis  functions  and  the  /“-closest  wavelet  basis  functions  are  not  "compactly  supported"  in 
the  intuitive  sense,  compared  to  the  Daubechies’  D4  wavelet  basis. 

Before  we  continue,  we  need  to  set  some  notations  used  in  the  rest  of  this  report.  Recall  our 
fixed  parameters.  Our  window  size  is  denoted  by  N,  where  N  =  2L  for  some  positive  integer  L. 
For  1  <  p  <  L,  we  will  abbreviate  v/:>'  by  /Vp.  Also,  we  shall  denote  the  Hilbert  spaces 

/■(Z.Vj.)  by  ‘Fp.  This  fits  with  our  previous  notation  of  *Fn  for  the  Hilbert  space  /“( Z,v).  The 
dimension  of  'Fp  is  Np.  We  will  use  the  words,  "vector."  "signal."  and  "function"  synonymously, 
although  "signal"  alludes  primarily  to  the  vector  being  written  in  the  standard  basis. 

Definition  1:  Let  f'i(n),  fi(n) . fp(n),  and  gp(n)  e  %  be  p+1  pair-wise  orthonormal 

vectors,  expressed  in  the  time  basis.  Consider  the  rotations  of  the  above  functions, 

(r~  kfm)(n)  =  fm(n  -  2mk),  for  k  =  0,  1 . Nm  -  1,  and  m  =  1,  2 . p,  including  the  rotations  of 

gp(n)  by  2pk,  k  =  0.  1 . Np~  1.  These  rotations  of  the  p+1  functions  make  up  /V  functions  in 

all.  We  say  this  set  of  N  functions  is  an  orthonormal  p-level  wavelet  basis  if  these  N  functions 
are  mutually  orthonormal  (and  hence,  span  an  /V-dimensional  space).  This  wavelet  baHs  is  also 
called  a  p+ Ist  generation  wavelet  basis.  Note  that  there  is  no  "scaling"  function  in  the  definition. 
See  References  (d [c|,  and  [k|. 

The  requirement  that  the  rotations  of  the  p+1  functions  be  orthonormal  is  written  as. 

(R;">f,n.Rr'kf„)  =  5,„.n  5,t. 

where  5X  x  is  the  Kronecker  delta  function.  One  may  also  replace  gp  for  fp  in  the  above  formula. 
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ELEMENTS  OF  THE  D4  WAVELET  BASIS 
CENTERED  WAVELET  FUNCTIONS  IN  TIME  SERIES 


So  far,  we  have  not  described  how  to  construct  such  a  wavelet  basis,  or  even  if  such  a  basis 
exists.  To  do  so,  we  begin  with  some  definitions.  Recall  the  notation  to  denote  a  vector  in  the  time 
basis:  f  =  fin).  We  may  also  use  the  boldface  (i.e.,  basis-free)  notation  of  the  vector  to  denote  the 
vector  in  the  time  basis. 


Definition  2:  A  vector,  fin]  €  ‘VL  _  j,  is  said  to  be  p-admissible  if 


(f(n  )|~  +  |f(n  +  /Vp)j"  =  2  for  all  n  e  Z^p. 


Definition  3:  A  pair  of  functions,  fin]  and  g[n]  <=  1;p  _i,  is  p-satisfactory  if  the  matrix 


J 


fin)  g(n) 

\  fin  +  Np)  g(n  +  Np) 


is  unitary  for  all  n  e  Note  that  each  function  of  a  p-satisfactory  pair,  (fin]  ,  g[n]),  is  p- 
admissible.  Compare  this  definition  of  p-satisfactory  to  Equation  (66)  in  Reference  [g]  and 
Equation  3.22  in  Reference  [fi. 

The  following  theorem  gives  a  constructive  method  for  making  an  orthonormal  wavelet 

basis. 


Theorem  1:  Suppose  yi,  <(>i  €  Tj).  Fix  a  p.  1  <  p  <  L,  and  define  yk,  0k  <=  H  - 1  by 
yk  =  Dk  l  yi  and  <jy  =  Dk  1  (j)  for  k  =  2,  3, ...,  p.  Let  =  0],  and  for  k  =  2,  3 . p  let 

fk  =  y,*Uy2*  ...  *Uk‘2yk.1*Uk'1<t>k 
and 

gr  =  yi*Uy2*  ...  *Ur'2yp.i*Up',yp. 


Then  the  collection  of  functions 


■  .  {R2Fkfp!0o  1  .  {R2r  kgp)k=o 


A/p  -  1 


is  an  orthonormal  p-level  wavelet  basis  for  %  if  and  only  if  yj,  0i  e  %  form  a  1 -satisfactory 
pair. 


2! 


The  rotations  of  gp  are  called  the  "extra"  -level  wavelet  basis  elements.  It  is  also  the 
p+lst  generaiion  of  the  wavelet  basis.  The  functions,  \|/t,  <J> x  e  Vo,  are  called  the  mother  and 
father  wavelet  functions. 

Example  1:  Haar  wavelet  basis. 


Let  0 1 (n ) 


VI 

‘VI 

o 


for  n  =  0 
for  n  =  1 
for  n  >  1 


VI 

and  \|/[(n)  =  J 

0 


for  n  =  0 
for  n  =  1 . 
for  n  >  1 


Then,  (<t>i,\j/i)  is  a  1 -satisfactory  pair.  Following  the  construction  given  in  the  above  theorem,  we 
get  functions  fj ,  f2-  ...,  fp,  gp  e  %,  which,  together  with  their  appropriate  rotations,  form  an 
orthonormal  basis,  called  the  p-level  Haar  basis.  By  following  the  computations  used  in  Theorem 
1 ,  one  may  check  that 

fi  =  .5  (1,  1,-1,  -1,0,0,  0 . 0), 

fi  =  j=(l,  1,  1,  1,-1, -1,-1, -1,0,0, ..., 0), 


and  so  on. 

Example  2:  Real  Shannon  Basis. 

Define  the  mother  and  father  wavelet  functions,  (jo  and  \j/i  e  V0  by 

if  0  <  m  <  A^2  -  1  or  3/V2  +  1  <  m  <  N  -  1 
if  m  =  /V2 
if  m  =  3/V2 

if  N2  +  1  <  m  <  3/V2  -  1 


VI 

~  i 

¥  i(m)  = 


and 
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<t>i(m) 


0  if  0  <  m  <  AS  -  1  or  3AS  +  1  <  m  <  N  -  1 
1  if  m  =  AS  or  m  =  3A/2 
il  if  AS  +  1  <  m  <  3N2  -  1  , 


so  that,  in  the  time  basis,  they  are  real-valued  and  the  pair  is  1 -satisfactory.  Follow  the 

A  /N 

construction  in  Theorem  1  to  get  the  wavelet  functions,  fj,  f2,  f3, ...,  fp,  gp,  where  f t  =  <J>  j  and 


and 


fk(m)  = 


2t/2  for  A/k  +  !  +  1  <  m  <  A/k  -  1 

orA/-Nk+l<m<Af-iVk+1-l 
2(k  -  l  )/2  for  m  =  Nk  +  i  or  m  =  N  -  /Vk  +  t 

i2(k- i)/2  form  =  ^k 

-  i2(k  -  1)/2  for  m  =  /V  -  /Vk 
0  otherwise 


,  k  =  2,  3 . p 


gp(m)  = 


2p/2  for  0  <  m  <  A^p  +  i  +  1 

or/V-/Vp+i  +  l<m</V  -  1 
;  j2(p-l)/2  form  =  /vp+i 

-i2(p ' 1)12  for  m  =  N  -  Np  +  i 
0  otherwise  . 


Note  that  the  supports  of  the  various  level  wavelet  functions  in  the  Fourier  basis  do  not  intersect  at 
more  than  two  points,  thus  displaying  the  strong  frequency  localization  of  the  real  Shannon 
wavelets.  (The  original  construction  created  a  high-pass  function  w  and  a  low-pass  function  v. 
Our  construction  creates  high-pass  and  low-pass  functions  for  different  levels.)  The  real  Shannon 
wavelets  are  also  spatially  localized.  See  References  [c]  and  [1]  for  more  details  of  the  real 
Shannon  wavelet  basis.  See  Figure  III- 1 . 


Example  3:  Daubechies1  D4  Wavelet  Basis 


This  wavelet  b  '  >s  generated  by  a  different  method,  compared  to  the  above  two 
examples  (and  the  /  - c., .  wavelet  basis).  We  do  not  actually  construct  the  wavelet  basis; 
instead,  we  create  the  wavelet  transform  so  that  we  can  take  the  wavelet  transform  of  any  signal 
function.  This  process  involves  constructing  two  sparse  matrices,  denoted  by  L  and  H  (for  Low 
and  High,  respectively).  The  sparseness  of  the  matrices  makes  the  computation  very  quick.  This 
is  why  the  algorithm  has  a  computation  time  of  O(N),  even  quicker  than  the  Fast  Fourier 
Transform. 
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The  matrices  are  made  up  of  four  different  numbers,  co,  cj,  ci,  and  c*.  These  are  exactly 
the  coefficients,  /i,v  («),  as  given  in  Table  1  of  Reference  [g],  for  N  =  2  (the  N  in  the  table).  Table  1 
of  Reference  [g]  only  gives  an  approximate  value.  For  N  <7,  these  coefficients,  h^in),  are 
computable  (they  are  algebraic  numbers).  For  our  case,  they  are 


/ 


(i  +  VT) 

4 

(3  +  11) 

4 

(3-V3) 

4 

(1  -  V~3) 

4 

0 


if  i  =  0 

if  i  =  1 

if  i  =  2 

if  i  =  3 
if  i  >  3 


We  then  create  the  matrices  from  these  coefficients.  Define 


U  = 


^2  i  -  j 


1  <=  i  <A^k 
1  <  j  i  N*  _ 1 


and 


Hk  = 


(— i  y + 'cj + 1  _  2  i 


1  <  i 

1  <j<iVk_, 


Also  define  Lk  =  Transpose  (Lk)  and  Hk  =  Transpose  (Hk).  These  matrices  are  defined  for  k  =  1 , 
2, ....  p,  where  p  is  the  level  of  the  wavelet  basis.  Now,  let  aP  =  f,  the  signal  (in  time  basis).  For  m 
=  p,  p  -  1, ....  2,  1,  define 


Then, 


am- 1  =  Lp_m ^  ia‘"  and  bm  =  Hp_m  +  iam. 


bp  1  contains  the  first  —  wavelet  coefficients, 

2 


bp  2  contains  the  next  ^  wavelet  coefficients. 


1 
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b°  contains  the  next  —  wavelet  coefficients,  and 

IP 

a0  contains  the  last  wavelet  coefficients. 

2P 

To  reconstruct  the  signal,  iterate  backwards.  Set 

am  =  Lp_m  +  jam  -  1  +  Hp_m  +  jbm  "  ^ 

for  m  =  1,2 . p.  Then,  af  =  f  as  desired. 

The  filter  bank  design  for  this  wavelet  transform  is  shown  in  Figure  III-4.  It  is  a  simpler 
design  than  the  other  wavelet  filter  banks,  since  it  only  involves  matrix  multiplication  with  a  sparse 
matrix.  The  sparseness  allows  for  a  fast  matrix  multiplying  algorithm. 

1.2.  Constructing  Wavelet  Bases 

The  construction  of  a  wavelet  basis,  as  given  in  the  above  theorem,  has  many  iterative 
steps.  This  suggests  a  possible  filter  bank  design  to  implement  this  wavelet  transform.  Indeed, 
we  have  developed  two  versions;  one  runs  in  the  time  domain  (see  Figure  III-5),  and  the  other 
runs  in  the  Fourier  domain  (see  Figure  III-6).  The  Fourier  domain  version  runs  substantially 
faster  than  the  time  domain  version,  due  mainly  to  fewer  convolutions  being  needed.  The  figures 
depict  a  2-level  wavelet  basis  transform. 
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Figure  III-4 

Daubechies'  D4  Wavelet  Filter  Bank  Design 


first  N/2 
coefficients 

next  N/4 
coefficients 

next  N/  2P 
coefficients 


last  N/  2P 
coefficients 


f 


=  Tf,  matrix  operation  on  a  vector. 
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Figure  III-6 

Fourier  Domain  Wavelet  Filter  Bank  Design 


Av,(n)  =  v,(-n)  n  =  0,  1 . N-  1 

v,  =  Fvj .  i 


\idft7 

~ \  / 


g 


=  Inverse  DFT  of  g 


There  are  other  ways  to  construct  wavelet  bases  (although  the  accompanying  filter  banks 
are  still  the  same  as  above).  The  following  theorem  gives  another  construction  process.  We  start 
with  a  definition: 

Definition  4:  For  a  function,  f  e  %,  define  Epf  €  Vp  for  p  =  1,2, ...,  L  by  the  following: 

2P-  1  ~ 

|f(n  +  k/Vp)|2,n  =  0,  U...,Np-\. 

~  k  =  0 

/v 

We  also  define  for  p  =  0:  Sof  =  f.  Note  the  definition  takes  a  function  written  with  respect  to  the 
time  basis,  uses  the  Fourier  transform  of  this  function,  and  yields  a  function  written  in  the  Fourier 
basis.  It  is  important  to  use  the  non-normalized  version  of  the  Discrete  Fourier  Transform.  This 
vector,  Lpf,  can  also  be  considered  as  a  vector  in  Vo  by  periodic  extension. 

Theorem  2:  Fix  a  level  p,  1  <  p  <  L.  The  following  are  equivalent: 


1.  (R2Pkf,f)  =  80.io  k  =  0,  1, ...,  Np  -  1.  Hence  (R2Pkf  )kP=  q  is  an  orthonormal  set. 

2.  (£pf)(n)  =  1,  for  all  n  =  0,  1 . /Vp-1. 

3.  f  =  V!  *  Uv2  *  U2V3  *  ■  ■  •  *  Up  "*vp  for  some  vj,  \2,  V3, ...,  vp  such  that  vk  is  k- 
admissible  for  each  k  =  1,2, ...,  p. 

4.  There  exists  a  p-level  wavelet  basis  whose  p^-generation  wavelet  functions  are 
rotations  of  f. 

Remarks: 

•  Statement  2  is  the  DFT  of  Statement  1 ,  hence  their  equivalence. 


The  vk's  of  Statement  3  are  determined  in  the  Fourier  domain  by: 


vk(n)  =  for  each  n  =  0,  1 . Nk  -  1 . 

(Ikf)(n) 
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•  To  create  a  p+1  generation  wavelet  basis  whose  p+l  st  generation  wavelet  functions  are 

rotations  of  f,  we  construct  wk,  k  =  1,2 . p,  such  that  (vk,wk)  is  a  k-satisfactory  pair 

for  each  k.  We  present  an  algorithm  for  this  process  later.  Now,  follow  the  procedure 
from  Theorem  1.  At  the  end,  we  get  fp  and  gp  with  f  =  gp.  To  get  f  =  fp,  just 
interchange  the  roles  of  fp  and  gp.  This  is  possible  since  both  functions  are  elements  of 

%■ 

•  Statement  3  can  also  be  written  in  the  Fourier  basis: 

f  =  Vj  •  Svi  •  S2V3  •  •  •  -  •  Sp ' 1  vp, 
where  ( •  )  denotes  the  outer  product. 

The  proof  of  Theorem  2  can  be  found  in  Reference  [m], 

1.3.  Constructing  k-Satisfactory  Pairs 

Let  v  e  _  i  be  a  k-admissible  function.  For  our  purposes,  we  assume  v  is  real-valued. 
We  would  like  to  find  a  w€  H-i  (also  real-valued)  so  that  (v,w)  is  a  k-satisfactory  pair. 
Assume  the  Fourier  transform  of  v  is  given  by: 

v(n)  =  |v(n)|  e'P<n>  (  By  convention,  take  p(n)  =  0  when  v(n)  =  0.) 

Then,  set  the  Fourier  transform  of  w  to  have  coordinates: 

w(n)  =  |v(n  +  (Vk)|  e'°(n). 


That  v  is  real-valued  implies: 
and 


p(0)  =  0  or  7t,  p(/Vk)  =  0  or  n 
p(n)  =  -  p(/Vk_  i  -  n)  (mod  2k) 


Similarly,  w  is  supposed  to  be  real-valued,  so  o(n)  satisfies  the  same  conditions.  The  condition 
for  a  pair  (v.w)  to  be  k-satisfactory  requires  an  accompanying  matrix  to  be  unitary.  By  assuming 
the  above  expressions  for  the  Fourier  transforms  of  v  and  w,  we  need: 

p(n)  -  p(n  +  /Vk)  =  k  +  CT(n)  -  o(n  +  /Vk)  (mod  27t)  for  n  =  0,  1 . Nk_  \. 
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The  choice  we  took  was  to  set: 


.  G(0)  =  p(0) 

.  a(Nk)  =  p(Nk)  +  n 

.  ct(n)  =  p(n)  -  k/2  for  n  =  0,  1,  2 . Nk  -  1 

*  o(n)  =  p(n)  +  n/2  for  n  =  Nk,  Nk  +  1, ...,  Nk_  i  -  1, 

all  modulo  2tl  This  is  the  procedure  we  used  to  generate  our  discrete  orthonormal  wavelet  bases. 

The  usefulness  of  wavelet  transforms  is  enhanced,  certainly,  by  their  multitude.  A 
particular  wavelet  transform  can  be  tailored  so  as  to  effect  a  desirable  response.  This  is  one  main 
advantage  of  the  discrete  wavelet  transform  over  the  Fourier  transform. 

2.  Closest  p-level  Wavelet  Bases 

The  variability  of  wavelet  bases  allows  one  to  choose  a  wavelet  basis  such  that  the  p^-level 
wavelet  function  fp  is  close  to  a  given  signal  s  e  %.  Specifically, 

Definition  5:  By  Statement  2  of  Theorem  2,  we  say  a  function  f  e  %  is  a  p^-level 
wavelet  function  if(Ipf)(n)  =  1  for  all  n  =  0,  1,  ...,  Np-  1.  We  can  also  extend  the  function 
(Zpf)(n)  periodically  to  be  a  function  on  (0,  1, ...,  N  -  1 }. 

Theorem  3:  Let  f  €  %.  Then  there  exists  a  p^-level  wavelet  function  fp  e  %  closest  in 
the  /2-metric  to  f  among  all  other  p^-level  wavelet  functions.  Moreover,  fp  is  given  by: 


fP(n)  = 


f(n) 

(£pf)(n)  ' 


Proof:  The  set  of  all  p^-level  wavelet  functions  is  a  compact  set  in  C^;  hence,  existence  is  clear. 

/V 

Now,  check  that  every  element  of  this  set  has  norm  equal  to  1.  The  function  fp  is  in  this  set  (one 
can  also  check  that  Xpf  is  a  p^-level  wavelet  function,  after  extending  periodically  to  be  a  function 

A 

on  {0,  1, ....  N  -  1 })  and  has  the  same  direction  as  f;  hence,  it  is  closest,  by  the  triangle  inequality. 

This  division  by  Zpf  is  what  we  shall  call  p-normalization.  For  p  =  L,  it  is  the  standard 
normalization. 
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The  closest  wavelet  need  not  be  unique.  Note  that  for  some  n  e  ((),  1,  ....  N  -  1 }  (Lpf)(n) 
may  be  zero.  In  this  case  then,  f(n)  must  also  be  zero;  hence,  fp(n)  is  undefined.  One  may  assign 

A 

fp(n)  to  take  on  any  value  so  long  as  fp  is  still  p-admissible.  The  option  we  took  when 

implementing  this  was  to  set  fp(n)  =  1  in  this  case.  This  allowed  us  to  generate  /2-ciosest  pth-level 
wavelet  functions  (and  hence,  p^-level  orthonormal  wavelet  bases)  to  any  function  f  e  %. 

A  particular  apolication  of  /“-closest  wavelets  construction  is  a  wavelet  basis  w'here  all  the 
wavelet  functions  from  each  level  are  close  to  a  given  function.  The  impetus  is  still  to  find 
orthonormal  bases  yielding  small  compression  loss  for  reproduction  and  detection  purposes. 

Theorem  4:  Suppose  f  e  %  satisfies  (Lpf)(n)  =  1.  Then  (Ip+if)(n)  =  1. 

The  proof  of  this  fact  follows  from  direct  computation  and  using  the  definition  of  the  expressions. 
Using  Theorem  4,  we  can  create  an  extremely  quick  algorithm  for  generating  a  p-level 
orthonormal  wavelet  basis  where  the  p^-level  wavelet  function  is  closest  to  a  given  signal. 

A  f  { n ) 

Let  fe  %  be  a  real-valued  signal.  Set  fp(n)  = - -,  w'here  Iif(n)  has  been  extended 

Iif(n) 

periodically.  Then,  v^n)  =fvi(n)|e1P(n' =  fp(n)  and  wpn)  =  \q(n  +  ^)e'°(n),  where,  as  before, 

p(n)  -  p(n  +  =  jc  +  a(n)  -  <T(n  +  ^-)  (mod  2^)  for  n  =  0.  1 . N  -  1.  Note  that  vfc  =  1  for  k  = 

2,  3 . p.  So,  creating  the  other  level  wavelet  functions,  f|.  fi . fp,  and  gp  is  much  easier; 


As  in  Theorem  2,  we  switch  the  definitions  of  fp  and  gp  to  get  the  desired  result.  The  only 
place  where  choice  is  involved  is  in  defining  the  w^'s.  We  used  a  very  simple  formulation.  For  k 
=  2,  3 . p,  define  wk  by: 

l  1  for  n  =  0 

for  0  <  n  <  Nk 
for  n  =  /Vk 
for  /Vk  <  n  <  /Vk  _  i  . 

We  now  give  another  construction  for  a  1 -level  orthonormal  wavelet  basis.  This 
construction's  main  purpose  is  to  correct  for  a  deficiency  in  using  the  standard  /  -closest  wavelet 
basis  for  compressed  detection  methods  (see  Chapter  IV-2).  Our  results  (see  Chapter  VII)  show 
that  the  /2-closest  wavelet  basis  performs  well  in  detecting  the  signal,  via  the  compressed  Neyman- 
Pearson  detector.  However,  if  the  signal  is  translated  by  one,  this  same  /"-closest  wavelet  basis 
does  not  detect  the  translate  as  well. 


Wk(n) 


-  i 


+i 


We  choose  a  1 -satisfactory  pair,  (f,g),  such  that  f  is  /2-closest  to  a  given  signal  \\i  and  g  is 
/2-closest  -  among  all  such  g’s  for  which  (f,g)  is  1 -satisfactory  -  to  Ry.  This  wavelet  basis 
induces  a  decomposition  algorithm,  called  the  ^-binary  decomposition.  Before  we  describe  the 
decomposition  algorithm,  we  first  derive  the  pair  (f,g)  for  vj t. 


Let  be  a  function  (in  the  time  basis)  in  %  =  /2( Z/v).  By  Theorem  3,  we  can  find  an  /"- 
closest  first-level  wavelet  function  to  \\r.  that  is,  there  exists  a  1 -admissible  function,  f,  which  is 
/2-c losest  to  V)/  among  all  1 -admissible  functions.  This  1 -admissible  function  is  then  the  mother 
wavelet  f,  and  the  corresponding  father  wavelet  g  is  determined  as  follows. 

Set  f(n)  =  ^-n)  -  ,  n  -0.  1 . N  -  1,  as  in  Theorem  3.  One  can  check  that  f  is  a  1- 

(IiVkn) 

admissible  function.  Moreover,  let  \j/(m)  =  a(rn)eia(m).  Then 


f(m)  = 


fla(m) 


V  a(m)2  +  a(m  +  N\)2 


giu(m) 
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in  the  standard  complex  polar  representation.  Note  that  the  phases  of  y  and  f  are  the  same.  As 
before,  we  construct  g  so  that  (f,g)  is  1 -satisfactory.  It  follows  that  the  Fourier  transform  of  the 
father  wavelet  g  has  the  form  g(n)  =  [f(n  +  N i  )|e'°^n)  and  that  a(n )  is  given  by  our  required  relation, 

G(m)  -  G(m  +  N\)  =  Jt  +  a(m)  -  a(m  +  N\)  (mod  2 tc),  for  m  =  0,  1,  ...,  N  -  1. 

Let  to  =  e10,  where  0  =  2k/N.  We  have  seen  that  (Ry)(m)  =  y(m)  wm,  where  R  is  the 
rotation  operator.  We  want  to  make  g  as  close  to  Ry  as  possible  in  the  /2-metric.  Therefore,  we 
want  to  minimize  the  value  D  given  by 


N-  1 

D  =  II  g  -  Ry  II2  =  II  g  -  Ry  If  =  ^  |  g(m)  -  y(m)  com  [ 2 

m  =  0 

/V,  -  1 

=  X  ( I  g(m)  -  y(m)  com  |2+  |  g(m  +  N\)  -  y(m  +  N\)  0)(m  +  iV|)  | 2) 

m  =  0 


N\  -  1 

s 

m  =  0 


VTa(m  +  N\)  e‘a(m) 
V  a(m)2  +  a(m  +  N\)2 

i2a.(m)  el0(m  +A,0 
Va(m)2  +  a(m  +  N  i)2 


-  a(m)  ei(a(ml  +  m0) 


-  a(m  +  N i)  e'(a(m  +  + 


2  \ 


Set  b(m)  = 


i2 

V  a(m)2  +  a(m  +  N\)2 


Thus,  we  need  only  minimize  each  term. 


A(m)  =  |  b(m)a(m  +  /V1)eio^m)  -  a(m)el(a(m)  +  m0)| “  + 

|  b(m)a(  m)el0,m  +;v')  -  a(m  +  Ni)e'<a(m  +/v'> +  m0  +  It> 


The  term  A(m)  is  minimized  when  the  arguments  of  the  terms  in  each  difference  are  identical,  i.e., 

o(m)  =  a(m)  +  m0 
and 

a(m  +  /Vj)  =  a(m  +  N\)  +  m0  +  k. 

Defined  in  this  way,  a  satisfies  the  phase  relation  above.  From  this  result  we  get. 
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g(m)  =  b(m)a(m  +  jVi)  e'(a(m'> +  mfl)  =  1 1  f(m  +N{)  ||  com. 


which  determines  g,  in  turn,  detennining  g.  Note  that  g  is  not  necessarily  the  /2-closest 
1 -admissible  function  to  Rv|/.  The  /2-closest  1 -admissible  function  to  R\y  is  Rf,  where 


Rf(n)  = 


\j/(n)  (on 
{IiHr)(n)' 


To  use  the  i)/-binary  decomposition  algorithm,  we  first  have  to  create  the  decomposition 
operator.  This  decomposition  operator  takes  a  vector,  s,  and  decomposes  it  into,  essentially, 
another  vector,  s  .  Note:  s  is  not  the  vector  s  written  with  respect  to  the  basis  generated  by  even 
translates  of  f  and  g. 


The  decomposition  operator  takes  the  vector  s  e  %  and  decomposes  it  into  two  vectors, 
Sf  and  sg,  lying  in  orthogonal  subspaces  of  equal  dimensions.  This  process  repeats  again  with 
these  two  vectors,  Sf  and  sg.  Specifically,  this  operator  projects  the  vector  s  into  two  orthogonal 

N 

subspaces  of  dimension  7,  respectively.  In  each  of  these  subspaces,  we  project  the  image  of  the 

original  vector  in  each  subspace  into  two  more  orthogonal  subspaces,  four  altogether,  and  each 
N 

having  dimension  4.  We  repeat  this  process  until  we  project  into  two-dimensional  subspaces. 
The  final  projection  into  one-dimensional  subspaces  is  canonical:  the  two-dimensional  subspaces 
are  spanned  by  a  particular  vector  (which  will  be  described  when  we  give  the  details  of  the 
operator)  and  its  orthogonal  complement;  thus,  the  projected  two-dimensional  vector  is  canonically 
projected  into  the  one-dimensional  subspace  spanned  by  that  particular  vector. 

Let  us  now  determine  the  first  set  of  subspaces.  Fix  a  real-valued  function,  \\i  (written  in 
the  time  domain)  in  %  =  /2( Z/v).  We  assume  this  function  is  normalized.  As  explained  above, 
we  get  the  72-closest  pair,"  f  and  g.  Both  f  and  g  are  also  real-valued.  The  even  translates  of  f  and 

g,  respectively,  span  mutually  orthogonal  subspaces  of  dimension  7.  Now,  project  the  function  \jr 
into  the  subspace  spanned  by  the  translates  of  f,  and  project  the  function  R\|/  into  the  subspace 
spanned  by  the  translates  of  g. 

Call  these  projected  functions  \|/f  and  \|/g,  respectively.  Now,  repeat  the  process  within 
each  subspace.  Find  the  "/2-closest  pair  to  \j/ f,"  fo  and  go,  and  the  "/2-closest  pair  to  \j/g," 
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f[  and  gi.  These  vectors  and  their  even  translates  span  four  subspaces,  each  of  dimension  j. 
Project  the  functions,  Vj/f,  R\|/f,  v| rg,  and  Rv|/g,  into  the  respective  subspaces  as  in  the  first  step 
above.  Repeat  the  process  until  the  subspaces  are  two-dimensional.  Note  that  this  operator,  once 
the  subspaces  are  determined,  is  a  linear  operator;  it  is  just  a  sequence  of  linear  projection  maps. 

This  set  of  nested  subspaces  forms  the  \|/-binary  decomposition  operator.  To  use  this 
operator  on  a  function,  s,  take  the  function  s  and  project  it  into  this  set  of  nested  subspaces.  The 
decomposition  of  the  projection  of  s  in  the  two-dimensional  subspaces  into  scalars  is  canonical;  the 
projections  into  the  two-dimensional  subspaces  give  pairs  of  numbers  that  will  be  the  scalars. 

This  projection  is  computed  taking  the  inner  products,  (s  ,  R2kf)  and  (s  ,  R2kg)  for 
k  =  0,  1 . —  -  1 .  These  inner  products  can  be  quickly  computed  by  convolution; 

N-  1 

(s,R2kf)  =  X  f(n-2k)s(n) 

n  =  0 

N-  1 

=  ^  f(n  -  2k)  s(n),  since  f  is  a  real  vector, 

n  =  0 

N-  1 

=  ^  (Af)(2k  -  n)  s(n),  where  (Af)(n)  =  f(— n),  the  involution  operator, 

n  =  0 

=  ((Af)*s)  (2k),  where  *  denotes  convolution, 

=  D((Af)*s)  (k),  where  D  is  the  downsampling  operator. 

The  inner  products,  (s  ,  R2kg),  are  computed  similarly.  Like  the  wavelet  filter  bank  designs,  these 
inner  products  can  be  computed  in  the  Fourier  domain  for  added  speed.  This  V(/-binary 
decomposition  is  illustrated  in  the  following  filter  bank  schematic  in  Figure  III-7.  Here,  the  first 
"/  -closest  pair”  will  be  denoted  by  f  and  g;  the  second  set  of  "/"-closest  pair"  will  be  denoted  by 
ff,  fg,  gf,  and  gg,  where  ff  and  fg  are  the  "/2-closest  pair"  in  the  subspace  spanned  by  f.  and  gf  and 
gg  are  the  "/2-closest  pair"  in  the  subspace  spanned  by  g. 
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y-Binary  Decomposition  Filter  Bank 


=  Dx.  (Dx)(n)  =  x(2n),  n  =  0,  1 . N  -  1 . 


CHAPTER  IV 


DETECTION  ALGORITHM 

In  this  chapter,  we  discuss  the  signal  detection  method  used  in  the  research.  The  method  is 
a  modification  of  the  Neyman-Pearson  (N-P)  detection  algorithm.  This  detection  algorithm  is 
optimal  in  a  certain  sense,  which  we  will  explain  later.  See  References  [n],  [o],  [1],  [m],  and  [c]. 

We  now  explain  the  N-P  detector  in  the  context  of  this  report;  that  is,  we  want  to  detect 
windowed,  discretized  signals.  Assume  we  receive  an  incoming  signal, 

x  =  (x(0),  x(l), ...,  x(N -  1)). 

We  would  like  to  determine  whether  this  incoming  signal  carries  the  prototype  signal  s  embedded 
in  noise,  where 

s  =  (s(0),  s(l) . s(/V-  1)), 

or  whether  the  incoming  signal  is  just  noise, 

n  =  (n(0),  n(l), ...,  n(N  -  1)). 


So  we  have  two  hypotheses; 


Hi  :  x  =  s  +  n. 


and 


Hq  :  x  =  n. 


Denote  the  choice  of  Ho  by  Co,  and  the  choice  of  Hi  by  Cj.  Then  there  are  four  possible 
outcomes, 

1.  We  correctly  identify  the  presence  of  the  signal.  The  probability  of  this  outcome  is 
called  the  probability  of  detection  (Pj),  denoted  by  P(C)  I  Hi). 
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2.  We  mistakenly  identify  the  presence  of  a  signal  when  in  fact  the  signal  is  not  present. 
The  probability  of  this  outcome  is  called  the  probability  of  false  alarm  (Pfa),  denoted  by 
P(C,  I  Ho). 

3.  We  incorrectly  determine  that  no  signal  is  present  when  in  fact  the  signal  is  present. 
The  probability  of  this  outcome  is  called  the  probability  of  a  false  negative  (Pn)  and  is 
denoted  by  P(Co  (Hi). 

4.  Finally,  we  correctly  determine  the  absence  of  the  signal  in  an  incoming  signal  that 
consists  of  just  noise.  The  probability  of  this  outcome  is  called  the  probability  of  a  true 
negative  (PIn),  denoted  by  P(Co  I  Ho). 

The  two  outcomes  with  which  we  are  most  concerned  are  the  first  two:  the  probabilities  of 
detection  and  false  alarm.  Moreover,  we  assume  an  active  sonar  system.  Thus,  we  should  have,  a 
priori,  the  probabilities  of  whether  the  signal  is  present,  P(Hi),  or  absent,  P(Ho).  These 
probabilities  will  affect  the  overall  outcome.  For  example,  there  is  an  expected  return  time  of  a 
transmitted  signal.  Around  this  expected  time,  the  probability  that  the  incoming  signal  carries  the 
prototype  signal  in  noise  increases.  Thus,  the  probabilities  of  detection  and  false  alarm  will  differ 
from  the  usual  situation  where  there  is  no  signal  expected. 

The  essence  of  the  N-P  detector  is  that  for  a  fixed  probability  of  false  alarm,  the  probability 
of  detection  is  maximized.  This  gives  the  optimality.  This  is  accomplished  through  computing  a 
maximum  likelihood  ratio.  For  a  threshold  (3  (determined  by  the  false  alarm  probability),  we 
choose  Hi  when  the  test  statistic  X,  is  greater  than  or  equal  to  the  threshold;  we  choose  Ho 
otherwise.  The  test  statistic.  X,  is  computed  by  the  maximum  likelihood  ratio: 

*  _pi(y) 

Po(y) 


where 


X  =  test  statistic, 

y  =  a  value  assigned  to  the  observation  x, 
p\(y)  =  probability  density  function  of  y  given  that  Hi  is  true, 
Pa(y)  =  probability  density  function  of  y  given  that  Ho  is  true. 
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The  probability  density  functions  are  Gaussian  normal  functions  with  means  and  deviations 
determined  empirically.  Essentially,  the  probability  of  false  alarm  sets  a  critical  value  for  y, 
deno'ed  by  yj.  This  critical  value  then  determines  the  value  of  the  threshold,  p. 

Without  going  into  details  (see  References  jn],  [o],  [m],  [c],  and  [1]  for  the  details  and  the 
computations),  the  test  statistic  is  computed  by  a  simpler  expression, 

X  =  sT  R-1  x, 

where  R  is  the  noise  correlation  matrix.  This  test  statistic  is  shown  (see  Reference  [c])  to  be 
independent  of  the  basis  for  expressing  the  vectors  and  matrix.  (In  this  sense,  this  implies  that  no 
particular  basis  choice  will  affect  the  outcome  of  the  detector.  In  particular,  wavelet  basis  or  the 
Fourier  basis  will  perform  just  as  well  as  the  rime  basis,  with  regard  to  the  N-P  detector.)  This  test 
statistic  is  essentially  an  inner  product.  If  the  noise  is  white  noise,  then  the  noise  correlation  matrix 
can  be  taken  to  be  the  identity  matrix.  We  will  assume  this  for  the  rest  of  the  report. 

The  drawback  to  using  the  full  N-P  detector  is  the  computing  time.  A  compressed  N-P 
detector  will  speed  up  the  process  by  computing  a  partial  inner  product;  that  is,  we  first  choose  a 
basis  for  computational  purposes.  We  then  take  M  (a  positive  integer  much  smaller  than  N ) 
coordinates  of  the  vectors  in  this  basis  and  compute  the  inner  product  with  only  these  coordinates. 
Here,  the  choice  of  the  basis  is  crucially  important,  as  can  be  seen  in  our  results. 

In  short,  let  A  be  a  subset  of  {0,  1, ...,  N  -  1 },  of  cardinality  L  and  B  be  a  subset  of  (0,  1, 
...,  N  -  1 }  containing  A  and  having  cardinality  M>L.  Then,  it  is  the  case  that  the  compressed  test 
statistic  based  on  the  partial  inner  product  with  coordinates  from  A  will  be  smaller  than  the  test 
statistic  based  on  the  partial  inner  product  with  coordinates  from  B.  Hence,  the  probability  of 
detection  grows  as  the  testing  subset  increases  in  cardinality.  That  is, 

Xa  <  Xu- 

On  the  other  hand,  choosing  a  "wrong”  subset  A  may  cause  Xa  to  be  small  relative  to  X.  It  may 
even  be  the  case  that  #(A)  >  #(B),  but  Xa  ^  provided  BdA.  Here,  #(A)  denotes  the 
cardinality  of  A. 

The  main  contention  is  that  when  compressed  to  a  subset  Be  {(),  1 . N  -  1 }  with 

respect  to  a  particular  basis,  it  is  always  the  case  that  the  compressed  test  statistic.  \r,  is  less  than 
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or  equal  to  the  full  test  statistic,  X.  Thus,  it  may  happen  that  X  >  (3,  yet  Xb  <  (3.  This  would 
increase  the  probability  of  a  false  negative  outcome. 

Cc  ntinuing  with  the  analysis,  one  constructs  the  compressed  Z-parameter  (see  Reference 
[m]),  which  in  our  case  is  given  as 


where  Ac{0,  1, ..., (V-l)  is  the  compressing  subset,  and  s  =  (s(0),  s(l ),...,  s(/V  -  1))  is  the 
signal  written  with  respect  to  some  fixed  orthonormal  basis. 

This  Z-parameter  is  related  to  the  test  statistic,  X.  The  test  statistic  gives  the  probability  of 
detection,  Pd,  as  a  function  of  the  probability  of  false  alarm,  Pfa.  The  graph  of  this  function  is 
known  as  the  Receiver  Operating  Characteristic  curve,  or  the  ROC  curve.  The  compressed  Z- 
parameter  arises  in  the  expression  of  Pd  as  a  function  of  Pfa, 

Pd  =  Q(Q-1(Pfa)-4), 

so  that  increasing  the  Z  values  improves  the  compressed  detector  toward  optimality.  (Here, 
Q(x)  =  1  -  <I>(x),  where  <I>(x)  is  the  standard  normal  distribution  function,  often  called  the  phi- 
function.  See  Reference  [m].  It  is  sufficient  to  know  that  Q  is  a  decreasing  function.  As  x 
increases,  Q(x)  decreases;  hence,  as  Z^  increases,  Q” 1  (Pfa)  -  Zg  decreases,  and  Pd  increases.) 

Finally,  we  measure  this  Z-parameter  as  a  dB  loss.  Define  the  degradation  factor  for  using 
the  compressing  subset  A  with  respect  to  the  fixed  orthonormal  basis  by 

D*  =  20  log10  Z£. 


Because  the  signal  is  normalized,  we  are  actually  computing  the  logarithm  of  the  ratio  of  the 
compressed  Z-parameter,  Z$,  to  the  full  Z-parameter,  which  is  1.  By  using  this  degradation 
factor,  we  can  see  how  much  energy  is  "lost"  due  to  compression.  We  often  plot  the  degradation 
factor  as  a  function  of  the  compression  ratio. 
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Thus,  no  degradation  occurs  when  there  is  no  compression  and  in  this  case,  we  have  0  dB  loss. 
We  reserve  describing  the  results  to  Chapters  V  and  VII. 

The  degradation  loss,  ?  can  also  be  interpreted  as  the  change  in  the  SNR  at  the  output  of 
detector.  Since  ranges  between  0  and  1,  ranges  from  -00  to  0,  respectively.  Optimality 
corresponds  to  =  0,  that  is  0  change  in  the  SNR. 


CHAPTER  V 


SIGNAL  SELECTION 

In  order  to  test  the  utility  of  wavelet  based  detectors,  we  created  a  set  of  "standard"  signals 
common  to  sonar  processing.  We  will  describe  them  in  detail  in  Section  1.  In  Section  2,  we 
describe  the  properties  of  these  signals.  Here,  we  will  also  discuss  the  Doppler  Ambiguity 
Displays  and  use  them  to  illustrate  the  differing  characteristics  of  the  signals.  In  Section  3,  we 
discuss  time-frequency  distributions  and  how  they  can  be  used  to  understand  signals.  Although 
we  have  specified  the  parameters  for  our  signals  (e.g.,  signal  length  and  sampling  frequency), 
these  parameters  may  be  adjusted  for  any  frequency  range  desired.  Indeed,  these  parameters  are 
user  inputs  in  out  C  signal  software  package. 

1.  Signal  Types 

We  study  five  types  of  signals:  continuous  wave  (CW),  hyperbolic  frequency-modulated 
(HFM),  linear  frequency-modulated  (LFM),  pulse  (PUL),  and  pseudo-random  noise  (PRN). 
They  were  all  computer  generated  using  the  Mathematica  mathematics  software  package.  The 
routines  are  now  currently  rewritten  in  C  for  faster  implementation  and  better  access  to  other 
detection  software  packages  like  the  multi-channel  detection  system  (see  Chapter  V).  These 
waveforms  are  described  in  the  appendix  of  Reference  [n]  and  are  also  described  in  References  [p] 
and  [q]. 

These  basic  parameters  for  all  these  signals  are  the  same: 

•  Sampling  Frequency,  fs  =  1000  Hz. 

•  Window  size,  N  =  1024. 

•  Signals  are  ail  normalized  (in  the  /2-norm)  to  have  unit  norm. 

•  Some  signals  (e.g.,  the  CHIP2  signal  as  described  in  Chapter  VI)  were  created  by 
taking  parts  of  the  above  signals  and  concatenating  them  together. 

Definition  of  the  signals  (n  =  0,  1, ...,  1023  =  N  -  1). 

•  CW:  s(n)  =  cos  (2  Ttfn  //s),  where  /  is  the  desired  frequency  (in  Hz),  a  user  input. 
The  equation  given  here  is  not  normalized. 
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HFM:  s(n)  =  cos 


-2  k  N 


AIl 


log 


Nfi  n 

(/2-/l)/s  U 


,  where  the  frequency  range 


(h  -/l)/s 

of  the  signal  modulates  from  the  lower  value  of  /'j  (in  Hz)  to  the  upper  value  of  f2  (in 
Hz);  they  are  user  inputs.  Again,  the  equation  here  is  not  normalized.  The  frequency  - 
vs.  -  time  plot  of  this  signal  should  yield  a  hyperbola  in  the  form  of  y  =  1/(C  -x), 
where  C  is  the  asymptote.  Aliasing  may  occur  if  the  frequency  range  extends  beyond 
the  maximal  discrete  frequency.  Note  that  the  upper  and  lower  frequency  bounds 
cannot  be  equal. 


LFM:  s(n)  =  cos 


rr  (2  n  -  /V  )[4  /V/m  +/w  (2  n  -  A/ )] 
4(V/s 


where  fm  is  the  median 


frequency,  and  fw  is  the  frequency  sweep,  both  in  Hz..  At  the  midpoint  of  time,  the 
signal  has  the  median  frequency.  The  frequency  sweep  from  the  beginning  of  signal  to 
the  end  of  the  signal  is  fw.  Again,  both  the  median  and  the  sweep  frequencies  are  user 
inputs.  The  frequency  -  vs.  -  time  plot  of  this  signal  should  yield  a  line  given  by  the 
equation  y  -  mx  +  b.  If  the  frequency  sweep  extends  beyond  the  maximal  discrete 
frequency  (Nyquist  frequency),  there  will  be  aliasing  as  in  the  HFM  case. 


•  PRN:  The  pseudo-random  noise  is  constructed  by  first  taking  a  completely  (uniform) 
random  noise  sample  (we  used  the  Random  [  ]  function  in  Mathematica)  and  taking 
its  Fourier  transform.  The  user  inputs  are  a  lower  and  upper  frequency  bound.  The 
lower  frequency  is  greater  than  or  equal  to  the  minimal  discrete  frequency  and  the 
upper  frequency  is  less  than  or  equal  to  the  maximal  discrete  frequency.  In  the  Fourier 
domain,  the  lower  and  upper  frequencies  are  lower  and  upper  bounds.  We  then  take 
the  Fourier  transform  of  the  random  signal  and  zero  out  all  the  discrete  frequency 
values  outside  the  interval  defined  by  the  upper  and  lower  frequency  bounds. 
(Actually,  because  the  random  signal  is  real-valued,  there  will  be  two  disjoint  intervals.) 
We  set  s(0)  equal  to  its  given  value  instead  of  zeroing  it  out  because  this  should  be 
nonzero. 


•  PUL:  s(n)  =  sin  (7t  n /e  //s)  sin  (2  k  n/S|g  / fs ),  where  is  the  carrier  frequency,  and  fsig 
is  the  signal  frequency,  both  in  Hz.  The  carrier  frequency  will  typically  be  much 
smaller  than  the  signal  frequency  so  that  half  the  period  for  the  carrier  is  equal  to  the 
window  length  (1.024  seconds  in  our  case).  Again,  both  frequencies  are  user  inputs. 
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The  definitions  of  the  above  signals  can  be  seen  in  the  schematics  of  a  time-frequency  plot 
of  Figure  V-l.  We  will  provide  more  detail  about  time-frequency  distributions  in  Section  3. 
There,  we  have  actual  Fourier  spectrograms,  frequency-versus-time  plots,  of  the  above  signals. 

Figure  V-2  shows  the  time-series  plot  of  the  above  signals.  Note  that  the  time-series  plot 
of  HFM  and  LFM  look  quite  similar:  the  differences  are  not  easily  discerned.  In  the  time- 
frequency  distributions,  these  differences  are  easier  to  discern. 

In  the  following  chapter,  we  use  a  "chip"  signal.  A  chip  -  also  known  as  a  chirp  -  is  a 
signal  built  up  by  concatenating  a  sequence  of  the  above  type  signals. 


I 


I 
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Figure  V-I 

Frequency  vs.  Time  Plots 


frequency  frequency 


time 


Continuous  Wave  HFK 
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CONTINUOUS  WAVE,  60  Hz 


2.1.  Properties  of  these  Signal  Types 


The  sonar  parameters  listed  in  Chapter  1  dictate  that  sonar  systems  operate  at  low 
frequencies  and  use  broad-band  and  long-duration  signals,  see  Reference  [q],  These  parameters 
determine  the  type  of  information  one  can  get  from  the  sonar  signal,  which,  in  turn,  is 
quantitatively  dictated  by  the  Doppler  Ambiguity  display  for  the  signal. 

A  Doppler  Ambiguity  display  shows  how  different  a  signal  is  when  translated  in  time  or 
Doppler-shifted  in  frequency,  or  both.  The  ambiguity  function  for  a  signal  s(r)  is  the  magnitude- 
squared  of  the  uncertainty  function.  One  defines  the  uncertainty  function  (the  uncertainty  comes 
from  changes  in  both  frequency  and  time;  compare  to  the  Heisenberg  Uncertainty  Principle  of 
quantum  mechanics)  by 


¥(x,8) 


s[(  1  —  5)  t]  s(t  +  t)  dr. 


The  variables  x  and  8  correspond  to  variation  in  time  and  frequency  respectively.  The  uncertainty 
function  measures  how  close  a  time-shifted  version  of  the  signal,  s(r  +  x),  is  to  a  frequency-shifted 
version  of  the  signal,  s[(  1  -  5)  r].  We  take  the  complex  conjugate  of  the  frequency-shifted  version 
because  the  signal  could  conceivably  be  complex-valued.  In  the  case  when  8  and  x  are  both  zero, 
then  we  get  the  I2- norm  of  the  signal.  This  is  where  the  uncertainty  function  achieves  its 
maximum  value.  The  ambiguity  function  for  the  signal  s(r)  is  then  given  by 


ATs(x.5)  =  Jv(x,S)|2. 


One  way  to  use  the  Doppler  Ambiguity  plot  is  to  take  the  return  signal  (provided  one  has 
detected  it)  of  a  prototype  signal  and  compute  the  inner  product  of  it  with  time-  and  frequency- 
shifted  versions  of  the  prototype  signal.  The  maximum  value  of  these  inner  products  will  give  the 
user  an  idea  of  how  much  shifting  is  due  to  time  (i.e.,  from  range  effects)  or  frequency  (i.e.,  from 
relative  motion  effects). 

Figure  V-3  shows  the  level  curves  of  a  hypothetical  ambiguity  function. 
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Figure  V-3 

Level  Curves  on  an  Ambiguity  Surface 


Our  standard  signals  all  produce  their  own  characteristic  ambiguity  diagrams.  We  will 
discuss  each  of  the  example  signals  separately.  The  examples  consist  of  two  types  of  each  of  the 
signal  types.  They  are: 

•  CW:  We  made  two  continuous  waves,  one  is  at  60  Hz  and  the  other  is  at  120  Hz. 

•  HFM:  The  first  hyperbolic  frequency'modulated  signal  has  a  wide  band;  it  sweeps 
from  20  Hz  to  200  Hz.  The  second  hyperbolic  frequency-modulated  signal  is  narrow- 
banded,  sweeping  only  from  80  Hz  to  120  Hz. 

•  LFM:  The  linear  frequency-modulated  signals  are  similar  to  the  HFMs.  The  first  is 
wide-banded,  with  a  linear  sweep  from  20  Hz  to  200  Hz,  and  the  second  is  narrow- 
banded,  sweeping  linearly  from  80  Hz  to  120  Hz.  Notice  some  aliasing  in  the  wide¬ 
band  signal,  near  the  upper  frequency  end. 
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•  PRN:  Again,  the  two  pseudo-random  noise  signals  are  banded  like  the  two  frequency- 
modulated  signals.  The  first  pseudo-random  noise  signal  is  banded  between  20  Hz  and 
200  Hz.  The  second  is  banded  between  80  Hz  and  120  Hz.  Note  that  the  second  PRN 
looks  less  stochastic  than  the  first,  although  both  were  randomly  generated. 

•  PUL:  Here,  both  pulse  signals  have  a  signal  frequency  of  100  Hz.  The  difference  is 
that  the  first  is  more  transient;  its  support  is  only  0.512  seconds.  The  second  pulse 
signal  has  a  longer  support,  taking  up  the  whole  window  (1.024  seconds). 

2.2.  Continuous  Wave  Signals 

Figure  V-4  shows  the  time-series  plot  of  the  two  continuous  waves.  Figures  V-5  and  V-6 
show  the  Doppler  Ambiguity  plots  for  the  two  signals.  Notice  in  the  Doppler  Ambiguity  plots 
that  the  main  lobe  is  an  almost  vertical  ellipse.  Also,  this  lobe  repeats  itself  in  the  time  (horizontal) 
direction,  suggesting  that  Doppler  shifts  can  be  measured  by  the  width  of  the  elliptical  lobe  (any 
one  of  them).  The  thinner  the  width,  the  greater  the  Doppler  shift.  One  can  determine  the  direction 
of  the  Doppler  shift  by  determining  whether  the  frequency  has  been  shifted  up  or  down. 

The  elliptical  lobe  shape  is  more  pronounced  in  the  higher  frequency  signal,  CW  120  Hz; 
see  Figure  V-6.  This  occurs  because  the  amount  of  frequency  shifting  due  to  a  fixed  relative  speed 
is  proportional  to  the  square  root  of  the  frequency.  Thus,  the  lobes  of  CW  120  Hz  are  shorter  in 
the  vertical  direction.  Thus,  the  CW  120  Hz  signal  is  better  for  Doppler  detection  than  the  CW  60 
Hz  signal. 

The  periodicity  in  the  time  direction  shows  that  these  signals  do  not  perform  well  in 
determining  time  differences  (i.e.,  range  determination).  If  one  sends  out  a  CW  signal  and 
compares  the  return  signal  with  time  translates  of  the  prototype  signal,  one  cannot  tell  whether  or 
not  there  have  been  delays  in  the  transmitted  signal.  As  seen  in  the  Doppler  Ambiguity  plots,  this 
drawback  is  even  more  pronounced  in  the  higher  frequency  signal  (another  example  of  the  effects 
of  the  Uncertainty  Principle). 
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2.3.  Hyperbolic  Frequency*Modulated  Signals 


Figure  V-7  shows  the  time-series  plots  of  the  two  HFM  signals.  Notice  the  second  signal 
is  almost  like  a  continuous  wave.  This  effect  is  also  seen  in  the  Doppler  Ambiguity  plots.  Figures 
V-8  and  V-9.  The  plot  for  the  wide-band  HFM  signal  shows  one  main  lobe  (which  has  a  fairly 
constant  width)  and  two  weaker  side  lobes.  The  constant  width  of  the  main  lobe  means  that  a 
Doppler  shift  of  the  original  signal  behaves  similarly  to  the  unshifted  version.  There  would  be  no 
side  lobes  if  the  signal  were  infinite  in  time  (or  at  least  sufficiently  long  in  time).  In  that  case,  one 
cannot  measure  Doppler  shifts  at  all.  On  the  other  hand,  one  can  measure  time  shifts  by  noting 
how  far  one  is  from  the  main  lobe.  For  an  infinite  length  signal  (i.e.,  no  side  lobes),  one  can 
measure  time  shifts  with  even  better  accuracy. 

The  narrow-band  HFM  signal  has  a  vague  main  lobe  shape,  consisting  of  many  small 
ellipses  (characteristic  of  CW  signals).  The  vagueness  of  the  main  lobe  says  that  although  one  can 
use  this  to  determine  time  shifts,  it  is  not  very  accurate.  Similarly,  one  can  also  determine  Doppler 
shifts,  but  not  as  well  as  when  using  CW  signals. 
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2.4.  Linear  Frequency-Modulated  Signals 


It  seems  as  though  there  is  not  much  difference  in  the  time-series  plots  of  linear  frequency- 
modulated  signals  versus  HFM  signals.  See  Figure  V-10  for  two  time-series  plots  of  LFM 
signals.  There  is,  however,  a  great  difference  in  their  Doppler  Ambiguity  plots.  One  can  also  see, 
from  the  next  section,  how  their  Fourier  Spectrograms  differ.  Again,  the  narrow-band  LFM  signal 
is  quite  similar  to  the  narrow-band  HFM  signal.  This  is  only  reasonable,  since  there  is  not  much 
room  for  the  frequency  to  modulate  in  either  case. 

As  seen  in  the  Doppler  Ambiguity  plot  of  a  wide-band  LFM  signal.  Figure  V-l  1,  one  can 
determine  both  time  and  frequency  shifts.  The  drawback  is  that,  aside  from  the  peak  at  the  middle, 
the  height  of  the  plot  is  uniformly  low.  This  puts  a  constraint  on  the  precision  of  the  measurement. 

The  narrow-band  LFM  signal,  as  mentioned  above,  is  very  much  like  the  narrow-band 
HFM  signal.  Even  the  Fourier  Spectrogram  distributions  are  not  much  different.  As  mentioned  at 
the  head  of  Subsection  2.1,  the  useful  choice  for  signals  would  be  wide-band  signals.  The 
similarity  of  the  two  narrow-band  signals,  HFM  and  LFM,  shows  why  wide-band  signals  are 
more  useful. 
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2.5.  Pseudo-Random  Noise  Signals 


The  first  pseudo-random  noise  signal,  band-limited  to  between  20  Hz  and  200  Hz,  looks 
more  like  noise  than  the  second.  Note  that  we  call  them  signals,  for  although  they  were  pseudo- 
randomly  generated,  they  can  be  re-created  exactly.  (A  pseudo-random  number  generator  is  a 
fixed  sequence  of  random  numbers,  a  random  number  generator  produces  random  numbers  that 
cannot  be  re-created.  They  are  truly  stochastic.  A  pseudo-random  noise  signal  is  not  only 
generated  by  a  pseudo-random  number  generator,  but  is  also  band-limited.  Thus,  pseudo-random 
noise  may  actually  be  a  signal  in  the  sense  that  it  can  be  re-created.  Figure  V-13  shows  the  two 
PRN  signals. 

The  Doppler  Ambiguity  plot  of  the  wide-band  PRN  signal.  Figure  V-14,  is  very  peaked  at 
the  center  and  drops  off  in  all  directions.  This  says  that  such  signals  have  poor  Doppler  or  time- 
shift  measuring  capabilities.  Their  use  is  primarily  for  covert  pinging.  The  randomness  hides  the 
signal  in  the  ambient  noise  of  the  medium,  and  the  signal  peak  of  the  Doppler  Ambiguity  plot 
means  one  can  get  good  estimates  for  range. 

The  narrow-band  PRN  signal's  Doppler  Ambiguity  plot.  Figure  V-15,  is  also  peaked,  but 
the  drop-off  is  not  as  sudden  as  in  the  wide-band  PRN  signal  Doppler  Ambiguity  plot.  This 
Ambiguity  plot  also  shows  just  a  hint  of  continuous  wave-like  properties.  One  can  measure,  with 
more  accuracy,  both  range  and  direction  changes  using  this  narrow-band  PRN  signal.  However, 
like  the  LFM  signals,  the  drop-off  only  makes  it  more  difficult  to  detect  differences. 
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2.6.  Pulsed  Signals 


Our  time-series  plots  of  the  two  pulse  signals  are  shown  in  Figure  V-16.  The  first  is  more 
transient  than  the  second.  Because  these  signals  have  a  constant  signal  frequency  (of  100  Hz),  they 
should  perform  somewhat  like  the  CW  signals.  This  means  they  should  be  able  to  distinguish 
Doppler  shifts.  Also,  because  these  signals  are  transient,  they  should  also  yield  some  range 
information.  To  some  extent,  they  do,  but  their  Doppler  Ambiguity  plots.  Figures  V-17  and  V-18, 
do  not  really  bear  this  out. 

For  the  short  pulse,  note  that  the  center  lobe  of  the  its  Doppler  Ambiguity  plot  extends  all 
the  way  from  the  bottom  to  the  top  with  uniform  thickness.  This  suggest  a  delta  function-like 
property.  A  delta,  or  spike,  function  will  yield  no  information  on  Doppler  shifts  and  will  have  a 
vertical  main  lobe  that  is  also  uniformly  thick. 

The  long  pulse  has  some  CW  signal-like  properties;  the  vertical  lobes  are  elliptical  in  shape, 
so  that  Doppler  shifts  can  be  determined.  We  believe  that  even  our  long  pulse  is  too  short.  A 
longer  pulse  within  a  larger  window  would  yield  much  better  information  on  both  frequency  and 
time  shifts. 
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3. 


Time-Frequency  Distributions 


In  this  section  we  will  briefly  describe  the  notion  of  time-frequency  distributions. 
Although  we  will  discuss  various  types  of  such  distributions,  we  will  illustrate  only  the  Fourier 
Spectrogram. 

The  basic  idea  of  a  joint  time-frequency  distribution  is  to  describe  the  energy  density  or 
intensity  of  a  signal  simultaneously  in  time  and  frequency.  From  harmonic  analysis,  we  know  that 
the  instantaneous  energy  of  a  signal  s(r)  (at  time  to)  is  the  square  of  the  absolute  value, 

js(f0)|2  =  intensity  per  unit  time  at  time  to. 

By  taking  the  Fourier  transform  of  the  signal  s (r). 


J' 

1*1 


s((o)  =  -Lr|  s(r)  e'1U)'  dr. 


we  can  also  determine  the  intensity  of  the  signal  per  unit  frequency: 


|s(w)p  =  intensity  per  unit  frequency  at  frequency  to. 


The  normalization  in  the  Fourier  transform  means  that 


j"  |s(t)pdr  =  |  |s(to)|2dto  = 


=  total  energy  =  1 


where  the  integrals  are  over  the  domain  of  the  respective  functions  and  we  will  set  the  intensity  of 
the  signal  so  that  the  total  energy  is  equal  to  1.  What  we  would  like  to  construct  is  a  joint  function 
of  time  and  frequency  that  represents  the  energy  or  intensity  per  unit  time  and  frequency.  We  have 
such  a  distribution,  namely, 

F(/,oo)  =  intensity  at  time  t  and  frequency  to. 


More  importantly,  integrating  over  all  time  or  over  all  frequencies  will  yield  the  intensity  per  unit 
frequency  or  the  instantaneous  intensity  respectively: 
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The  total  energy  of  the  signal  is  then  expressed  as  the  integral, 

Energy(s)  =  j  F(?,co)  dco  dr. 

This  joint  function  of  time  and  frequency  would  be  our  distribution.  Its  plot  is  a  surface 
given  as  a  function  of  two  variables,  time  and  frequency.  Our  plots  are  given  in  density  plot 
format.  That  is,  the  darker  portions  correspond  to  stronger  intensities.  These  plots  for  the  various 
example  signals  are  given  in  their  respective  sections. 

There  are  a  number  of  time-frequency  distributions.  The  most  notable  ones  are  the 
Wigner-Ville  Distribution,  the  Rihaczek  Distribution,  the  Page  Distribution,  and  the  Fourier 
Spectrogram.  We  are  currently  developing  a  Wavelet  Spectrogram,  which  performs  better  than 
the  above  in  that  it  does  not  need  to  iterate  over  subwindows.  Our  Wavelet  Spectrogram  uses  the 
fact  that  the  wavelet  transform  already  possesses  time  and  frequency  information.  The  drawbacks, 
at  the  moment,  are  that  the  Wavelet  Spectrogram  does  not  yield  as  good  a  frequency  resolution  as 
one  would  like  and  that  the  time  resolution  is  only  relative  (i.e.,  one  needs  to  do  additional 
computations  to  determine  real-time  values).  One  advantage  of  the  Wavelet  Spectrogram  is  the 
computational  speed.  Another  advantage  is  that  the  choice  of  wavelet  basis  will  allow  for  better 
frequency  resolution. 

We  illustrate  the  time-frequency  distributions  with  the  Fourier  Spectrogram  (also  known  as 
the  windowed  Fourier  Spectrogram  or  the  time-dependent  Fourier  Spectrogram)  for  the  signals  we 
created  in  the  previous  section.  Note  the  difference  in  the  HFM  and  LFM  signals.  It  is  clear  in 
these  illustrations  that  the  frequency  is  linearly  modulated  for  the  LFM  signal.  It  is  not  so  clear  that 
the  frequency  is  hyperbolically  modulated  for  the  HFM  signal,  since  the  modulation  is  fuzzy  near 
the  end  (due  to  the  finiteness  of  the  signal). 


The  PRN  signal  is  seen  to  have  a  fairly  uniform  distribution  within  the  frequency  ranges 
set  by  the  band  limits.  This  indicates  the  random  number  generator  is  a  uniform  random  number 
generator  (as  opposed  to  a  normal  random  number  generator). 

Finally,  note  that  the  ends  for  the  PUL  signals  are  cleaner  than  the  ends  of  the  CW  signals. 
This  is  because  we  use  the  rectangular  window  function  to  create  the  subwindows  for  the  Fourier 
transforms.  Thus,  for  the  CW  signals,  we  get  an  abrupt  drop-off,  whereas  the  PUL  signals  yield 
gradual  drop-offs. 

To  generate  the  Fourier  Spectrograms,  we  set  a  subwindow  size  M  (in  our  case,  M  =  12X) 

and  a  step-size  r  (  =  X).  Let  s  =  s[n|  =(s(0),  s(l) . s(1023)  )  be  the  time-series  expression  of 

the  signal.  Then,  the  distribution,  F(/,co),  is  determined  by  a  Fourier  transform: 

M  - 1 

F (r,(o)  =  £  s(tr  +  m)  e~  2  711  m 

m  =  0 

where  t  =0,  1 . and  w  =  (),  1, ...,  A/  -  1.  Because  of  the  symmetry  in  the  frequency 

direction  (since  our  signals,  s,  are  real-valued  in  time),  we  take  only  the  first  half  of  the  discrete 

frequencies:  co  =  0,  1 . ^  .  This  is  why  the  vertical  frequency  axis  ranges  from  0  to  65. 
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Fourier  Spectragram  of  CW  60  Hz 
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Fourier  Spectragram  of  CW  120  Hz 


i)Z-A  ajnSj j 


Aouenbejj  ajajosiQ 
\Z-\  ajn3jj 


Discrete 


9  L 


/ouanbajj  ajajosig 
ZZ-\  3Jn3u 


Fourier  Spectragram  of  LFM  20-200  Flz 
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Fourier  Spectragram  of  LFM  80-120  Hz 
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Fourier  Spectragram  of  PRN  20-200  Hz 


Fourier  Spectragram  of  PRN  80-120  Hz 
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Fourier  Spectragram  of  PUL  100  Hz  (Short) 


CHAPTER  VI 


THE  SIMULATION  SYSTEM  AND  A  DESCRIPTION  OF  THE  EXPERIMENTS 

In  this  chapter,  we  describe  the  transmission  simulation  system  that  we  built  to  test  the 
efficacy  of  the  active  signals  we  created.  However,  the  primary  emphasis  of  this  chapter  is  to 
describe  the  experiments  using  the  transmission  simulation  system.  The  results  of  the 
experiments  and  the  data  interpretation  are  given  in  the  following  chapter.  Section  1  of  this  chapter 
will  explain  each  of  the  modules  that  composes  the  simulation  system.  Section  2  explains  our 
multi-channel  detection  algorithm,  which  aims  to  detect  families  of  signals,  classified  by  possible 
Doppler  shifts  or  translational  shifts.  Section  3  presents  an  example  of  the  compressed  Neyman- 
Pearson  detection  algorithm  (see  Chapter  IV)  implemented  on  a  hyperbolic  frequency-modulated 
signal.  CH1P2.  This  section  shows  how  this  signal,  CHIP2,  changes  due  to  Doppler  and  bottom- 
bounce  effects.  Moreover,  Section  3  also  shows  that  the  /2-closest  wavelet  basis  outperforms  all 
other  bases  with  regard  to  the  compressed  N-P  detector.  Section  4  will  describe  the  experiments 
using  the  multi-channel  detector  on  the  signal,  CHIP2. 

1.  The  Transmission  Simulation  System 

The  transmission  simulation  system  is  also  described  in  the  User’s  Guide  (Reference  [q]). 
The  User's  Guide  is  intended  as  a  tutorial,  and  thus,  it  does  not  give  much  motivation  or 
explanation  as  to  how  is  was  built. 

Overall,  the  transmission  simulation  system  (TSS)  consists  of  eight  modules.  Each 
module  is  a  separate  routine  that  simulates  a  particular  phase  in  the  transmission  history  of  the 
signal.  These  modules  are: 

•  Signal  Generator 

•  Projector  Response 

•  Transmission  (multi-path) 

•  Reflection 


Doppler 


Reverberation 


•  Hydrophone  Response 

•  Receiver/Detector 

The  first  module,  the  signal  generator  module,  creates  the  signal.  This  module  is  not  technically 
part  of  the  TSS.  The  detector  part  of  the  last  module  is  also  not  technically  part  of  the  TSS.  These 
two  are  actually  computed  outside  of  the  system.  In  practice,  as  currently  implemented,  we  have  a 
separate  routine  that  generates  the  prototype  signal  and  a  separate  routine  that  receives  signal-like 
data  and  runs  the  compressed  Neyman-Pearson  detector,  with  an  appropriate  transform  as  its  only 
user  input. 

1.1.  Signal  Generator 

The  signal  generator  module  creates  signals  like  those  described  in  Chapter  V.  Such 
signals  begin  with  norm  equal  to  one  (in  most  cases).  They  are  then  amplified  so  that  the  signal- 
to-noise-ration  (SNR)  is  about  230  dB.  (The  230  dB  is  a  default  value;  the  user  can  set  any  SNR 
value.)  Setting  the  SNR  level  is  part  of  the  TSS,  hence  we  put  this  amplification  phase  among  the 
simulation  modules.  The  amplification  routine  is  a  straightforward  linear  scaling.  In  addition  to 
scaling,  we  also  add  two  windows'  worth  of  noise. 

The  noise  is  actually  zero  noise,  meaning  that  we  take  the  signal  -  a  finite  sequence  of 
length  /V  =  1024  -  and  pad  this  sequence  in  front  and  back  with  2048  zeros.  The  output  of  the 
signal  generator  is  then  a  signal  of  length  5  N  =  5120.  The  reason  for  the  padding  is  that  the 
transmission  module  is  a  finite  impulse  response  filter  that  needs  at  least  two  windows'  worth  of 
data.  We  need  to  simulate  the  effects  of  the  transmission  also  on  noise,  in  addition  to  the  effects 
on  the  signal. 

1.2.  Projector  Response 

The  projector-response  module  simulates  the  nonflatness  of  most  projectors.  Projectors, 
because  of  their  physical  construction,  do  not  project  with  equal  intensity  over  all  frequencies.  The 
intensity  at  each  frequency  can  be  graphed  as  a  function  of  the  frequency.  Such  graphs  are  called 
frequency  response  curves. 
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Most  projectors  can  be  built  to  have  a  fairly  flat  response  curve  over  the  significant  range. 
However,  at  the  intensity  levels  of  230  dB  SNR.  there  will  be  significant  distortions,  and  the 
response  curve  will  not  be  flat.  Often,  the  upper  frequency  end  will  drop  off,  relative  to  the  whole 
frequency  spectrum.  Since  we  are  engaged  in  creating  low-frequency  active  signals  (i.e., 
frequency  range  between  20  Hz  to  500  Hz),  the  drop-off  at  the  high  end  is  still  small  even  at  230 
dB  SNR.  This  suggests  that  the  flat  response  default  (which  we  use)  is  a  suitable  default. 

The  input  of  this  module  is  the  scaled  signal  of  length  5  N  =  5120.  The  output  is  also  of 
length  5  N  =  5120.  The  routine  simply  multiplies  each  entry  of  the  sequence  by  a  scaling  factor 
corresponding  to  the  response  curve  and  returns  this  value  to  the  same  position  in  the  sequence. 
For  the  default  setting,  the  scaling  factor  is  set  at  1  for  all  entries. 

1.3.  Transmission  Effects 

The  next  phase  after  the  projector  response  is  the  effect  simulated  by  the  transmission 
module.  This  module  simulates  the  effect  of  bottom  bounce  encountered  in  a  shallow-water 
channel  scenario.  The  result  of  bottom  bounce  is  that  the  received  signal  is  a  conglomeration  of  the 
signal  taking  various  paths.  This  effect  is  also  known  as  multi-path. 

The  module  routine  is  a  finite  impulse  response  filter,  meaning  that  it  takes  a  finite 
sequence  and  returns  a  finite  sequence  (possibly  of  a  different  length).  We  will  first  describe  the 
model  upon  which  the  module  is  based.  Then,  we  will  give  the  equations  for  the  filter. 

This  module  employs  the  Pecris  model  for  bottom  bounce  effects,  as  illustrated  in  Figure 
VI- 1.  When  the  sonar  signal  impinges  on  the  bottom,  there  is  attenuation  loss  because  of  the 
change  in  the  signal  medium  from  water  to  solid.  The  amount  of  attenuation  loss  is  dependent  on 
the  angle  of  incidence,  a,  measured  from  the  surface,  as  well  as  the  specific  gravity  of  the  sea  bed 
and  the  water  pressure  at  the  bottom.  However,  there  is  a  critical  angle.  aL,  such  that  for  angles  of 
incidence  greater  than  ac.  the  signal  will  be  completely  absorbed  by  the  bottom.  This  sets  a 
maximum  number  of  bounces  off  the  bottom,  given  as  a  function  of  the  distance  between  the 
source  and  target,  their  respective  depths,  the  depth  of  the  water  channel,  and  the  critical  angle.  ac. 
There  is  no  attenuation  at  the  air-to-water  interface. 

The  existence  of  a  maximum  number  of  bounces  allows  us  to  model  this  as  a  Finite 
Impulse  Response  (FIR)  filter.  In  this  filter,  the  signal  (of  size  5  N)  is  entered,  and  the  output  is  a 


longer  duration  signal  (the  size  is  about  IN)  because  the  multi-reflected  paths  are  longer.  The 
earlier  arriving  signal  has  less  attenuation  than  the  later  arriving  signal.  The  amount  of  attenuation 
is  empirically  calculated  to  be  about  45  dB/sec.  when  given  the  parameters  for  the  density  of  the 
sea  bed  and  the  respective  speeds  of  sound  in  water  and  in  the  sea  bed.  The  Pecris  model  and  its 
equations  were  supplied  by  DARPA/AT&T  Acoustic  Warfare  research. 

The  equation  for  the  filter  is 

M  -  i 

s'(n)  =  £  s(n  —  j )  h(j). 

j  =  o 

where  h[n|  =  (h(0).  h(  1 ) . h(M  - \))  is  the  filter,  s|n]  =(s(0),  s(l) . s(5  N  - \))  is  the 

incoming  signal,  and  s'|n)  =(s’(0),  s'(l) . s'(3  N  -  1))  is  the  output  signal.  The  length  of  the 

filter,  M,  is  usually  2  N,  because  of  the  channel  depth  and  the  source-to-target  range.  This  is  why 
the  original  signal  was  padded  by  two  windows'  worth  of  zero  noise.  To  be  precise,  one  should 
index  the  output  signal,  s'|  n  |,  as  s'(2  N  -  1 ),  s'(2  N) . s'(5  N  -  1 ). 

Note  that  the  output  is  only  three  windows  long.  The  filter  does  not  shorten  the  signal. 
Indeed,  it  actually  lengthens  the  actual  signal,  which  has  length  A/,  to  3  N.  The  output  of  the  filter  is 
all  signal,  with  the  first  and  last  entry,  s'(0)  and  s'(3  N  -  1),  equalling  zero.  This  output  is  then 
padded  again  by  two  windows'  worth  of  zero  noise  in  both  the  front  and  the  back.  Thus,  the  final 
output  of  the  module  is  a  signal  of  length  IN. 


Figure  Vl-1 
Bottom  Bounce  Model 

about  20  nm 


The  model  should  yield  a  linear  drop  in  signal  strength  (measured  in  dB's)  as  a  function  of 
time.  Because  the  arrival  of  the  reflected  signal  occurs  at  discrete  moments  of  time,  the  FIR  filter 
is  finite  only  at  those  discrete  time  moments  corresponding  to  the  arrival  of  the  reflected  signal,  at 
all  other  points  it  is  minus  infinity.  (Note  that  minus  infinity  intensity  in  dB  is  equal  to  zero 
intensity  measured  as  amplitude  height.) 

This  module  is  used  again  after  the  Doppler  module.  This  is  to  model  the  returning  portion 
of  the  signal's  trip  from  source  to  target  to  receiver.  In  the  second  use  of  this  module,  the  input- 
signal  length  is  still  7/V.  The  output-signal  length  is  then  9 N.  This  consists  of  two  windows  of 
noise  followed  by  the  signal,  stretched  out  over  five  windows  due  to  this  bottom  bounce  effect, 
and  finally  followed  by  two  more  windows  of  noise.  The  noise  window  is  just  zero  noise,  so  the 
upshot  i.-,  five  windows  of  signal  padded  in  front  and  back  by  two  windows'  w  orth  of  zeros  each. 

1.4.  Reflection 

This  module  uses  either  the  monostatic  target  strength  or  the  average  cosine  approximation 
to  the  target  strength  for  reflection  of  signals  off  targets.  This  model  is  also  a  finite  impulse 
response  filter.  The  input-  and  output-signal  lengths  are  equal:  7 N.  See  Reference  [o]  or  |r]  for 
detaded  descriptions  of  the  reflection  model.  Our  default  input  is  the  monostatic  target  strength 
model. 

There  was  no  change  of  this  module  from  the  previous  presentation  of  the  TSS  (see 
Reference  [a|). 

1.5.  Doppler 

The  Doppler  module  is  also  a  finite  impulse  response  filter.  It  filters  via  a  sine  function 
(i.e.,  sim nit  function).  The  output  of  the  Doppler  module  may  spill  the  signal  outside  the 
window  .  This  occurs  when  there  is  a  negative  Doppler  shift  (i.e.,  the  target  is  moving  away  from 
the  signal  source/receiver).  Nevertheless,  the  length  of  the  output  signal  from  th:.,  module  is  equal 
to  the  length  of  the  input  signal,  w'hich  in  our  case  is  still  7/V. 

The  Doppler-shift  filter  is: 

L  - 1 

Sjj,  p(n)  =  V7  s(k)  SincOi  t  -  k). 

k  =  ii 

XX 


where  s[n]  is  the  incoming  signal  of  length  L ,  t  is  the  amount  of  Doppler  shift,  and  Sinc(.v)  is  the 

sine  function:  Sinc(.v  )  =  Sm^A  \  The  effects  of  Doppler-shifting  a  signal  are  shown  in  several  of 
the  time-series  plots  of  Doppler-shifted  versions  of  our  example  signal,  CHIP2. 

1.6.  Reverberation 

This  module  models  reverberation  noise  as  a  moving-average  stochastic  process.  At  any 
given  time,  reverberation  from  each  sample  of  the  signal  is  occurring  at  some  random  intensity. 

The  moving  average  equation  is  also  a  FIR  Filter: 

M  -  1 

v(m)  =  X  a(m  - 

k  =  0 

where  M  is  the  length  of  the  signal  (in  our  case,  M  =  9N  ),  s[k]  =  (s(0),  s(l) . s(M  -  1))  is  the 

signal  (with  the  zero  paddings  in  front  and  back),  and  a[k]  =(....  a(-l),  a(0),  a(l),  a(2), ...)  is  an 
infinite  sequence  of  normal  random  draws  with  mean  =  0.0  and  standard  deviation  o  =  input 
value.  The  sum,  v[mj,  is  then  the  noise.  This  noise  is  added  to  the  signal  and  the  resulting 
sequence  is  the  output  of  the  module. 

There  is  no  default  value  for  this  module. 

1.7.  Hydrophone  Response 

The  function  and  code  for  this  module  is  identical  to  the  projector-response  module. 

1.8.  Comments 

The  actual  filter  codes  used  in  the  TSS  are  based  on  the  equations  given  in  the  above 
subsections.  The  actual  algorithms  are  modifications  of  these  equations  due  to  computing  time 
and  overflow/underflow  constraints. 

The  following  chapter  includes  illustrations  of  the  time-series  plot  of  a  signal  (CH1P2) 
affected  by  bottom  bounce  and  Doppler  effects. 
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2.  The  Multi-Channel  Detector 

The  multi-channel  compressed  Neyman-Pearson  detector  is  a  system  designed  to  detect  the 

presence  of  any  translate  and  Doppler  shift  of  a  returned  signal.  We  begin  with  a  prototype  signal 

(called  "CHIP2"  in  our  experiments)  of  window  size  N  =  1024.  This  signal  is  zero-padded  at  one 

•) 

end  to  make  a  signal  of  size  2N.  This  padded  signal  is  used  to  create  a  p-level  /“-closest 
orthonormal  wavelet  basis.  This  is  done  using  the  construction  described  in  Theorem  III-3. 

To  create  the  multi-channel  detector,  we  first  generate  Doppler  shifts  of  the  original  signal, 
from  0  kts  to  ±m  kts  (m  stands  for  maximum).  Each  distinct  Doppler  shift,  including  the  zero- 
shifted  signal,  corresponds  to  a  channel  of  the  multi-channel  detector.  Each  Doppler-shifted  signal 
is  also  zero-padded,  similar  to  the  non-Doppler-shifted  signal.  We  have  a  separate  Doppler- 
shifting  code  that  produces  the  Doppler  shift  of  an  input  signal.  This  code  is  similar  to  the  code 
used  in  the  Doppler  module  of  the  TSS. 

Now,  each  channel  (i.e.,  each  distinct  zero-padded  Doppler  shift  of  the  signal)  is  translated 
by  one.  N  times.  This  translation  is  done  to  model  the  effect  of  not  capturing  the  returned  signal 
properly  in  the  receiving  window.  For  each  channel  and  each  translate,  we  compute  the  wavelet 
coefficients  with  respect  to  a  fixed  wavelet  basis.  The  wavelet  basis  is  now  the  detector  in  the 
sense  of  the  compressed  Neyman-Pearson  detector  (see  Chapter  IV). 

The  wavelet  basis  is  generated  by  the  wavelet  function  that  is  /“-closest  to  the  original 
signal.  In  some  of  our  experiments,  we  also  used  wavelet  bases  generated  by  wavelet  functions 
that  are  /2-closest  to  the  Doppler  shifts  of  the  original  signal,  as  well  as  multi-path  versions  of  the 
original 


Fix  a  compression  size  M  <  N.  This  number  is  the  size  of  the  subset  of  all  the  wavelet 
coefficients  that  we  will  use  in  the  compressed  Neyman-Pearson  detector.  For  this  fixed  M,  we 
choose  the  M  largest  weights  of  the  wavelet  coefficients  for  each  channel  and  each  translate.  This 
will  determine  the  subset  of  the  wavelet  coefficients  used  in  the  compressed  detector. 

Because  of  the  recursive  nature  of  the  wavelets,  all  these  data  are  actually  quite  easy  to 
store.  It  is  not  essential  to  store,  for  example,  all  translates  of  each  channel;  the  first  2P  translates 
(including  the  zero  translate)  are  sufficient  to  supply  the  data  for  the  N  translates. 
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We  now  receive  a  signal.  In  other  words,  we  send  a  prototype  signal  (e.g.,  a  CW  or  an 
LFM  signal)  through  the  TSS  and  take  the  output  of  the  simulator  as  the  return  signal.  This 
returned  signal  will  look  significantly  different  from  the  prototype  signal  that  we  send  out  (or 
rather,  into  the  simulation  system).  It  may  also  be  the  case  that  the  received  signal  is  just  noise;  the 
prototype  signal  is  not  present  in  that  window. 

Take  the  wavelet  transform  of  the  returned  signal,  using  the  same  wavelet  basis  as  the  one 
chosen  to  create  the  detector.  Use  the  subset  determined  above  to  check  for  the  presence  of  the 
prototype  signal  in  the  received  signal.  We  then  check  each  channel  to  determine  whether  Doppler 
shifting  has  occurred  in  the  signal. 

A  diagram  of  the  multi-channel  detector  methodology  is  given  in  Figure  VI-2. 
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Figure  VI-2 

The  Multi-Channel  Detector  System 
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When  using  the  multi-channel  detector  on  the  output  of  the  TSS.  the  detector  will  scan 
through  the  nine  window-lengths,  yielding  8  N  data  points.  (It  is  not  necessary  to  scan  the  very 
last  point  of  the  ninth  window,  so  there  are  only  8  whole  windows  to  consider.  )  A  typical  output  of 
the  multi-channel  detector  will  be  described  at  the  end  of  the  following  example  section. 

3.  The  Example  Signal 

We  created  a  signal  called  CHIP2.  This  signal  is  the  concatenation  of  five  HFM  signals. 
The  sampling  frequency  is  1000  Hz,  and  the  window  size  is  1024;  so  the  signal  has  a  duration  of 
1.024  seconds.  The  HFM  signals  are: 

•  The  first  0.27  sec.  is  HFM  ranging  from  80  Hz  to  120  Hz. 

•  The  second  0.24  see.  is  HFM  ranging  from  120  Hz  to  180  Hz. 

•  The  third  0.24  sec.  is  HFM  ranging  from  180  Hz  to  100  Hz. 

•  The  fourth  0. 10  sec.  is  HFM  ranging  from  100  Hz  to  150  Hz. 

•  The  fifth  and  last  0.17  sec.  is  HFM  ranging  from  150  Hz  to  80  Hz. 

See  Figure  VI-3  for  a  time-series  plot  of  CHIP2.  (Figure  VI-3  and  all  following  figures  are  given 
after  the  text.) 

•*) 

In  using  the  basic  compressed  Neyman-Pearson  detector,  we  create  the  /"-closest  3-level 
wavelet  basis  based  on  the  signal.  See  Figure  VI-4  for  the  time-series  plot  of  the  /2-closest  3r^- 
level  wavelet  function.  Figure  VI-5  shows  the  difference  between  the  CHIP2  signal  and  the 
/2-closest  3rd-level  wavelet  function.  Figure  VI-6  shows  the  time-series  plots  of  the  wavelet 
functions  at  all  levels:  first,  second,  third,  and  the  extra  third-level  wavelet  functions. 

We  also  sent  this  signal  through  our  transmission  simulation  system,  simulating  a  +5  kt 
Doppler,  a  shallow-channel  bottom  bounce  effect  (multi-path),  and  the  combination  of  both:  +5  kt 
Doppler  and  the  effect  of  multi-path.  Figure  VI-7  shows  the  time-series  plots  of  CHIP2  under 
these  scenarios,  compared  to  the  prototype. 
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4.  The  Experiments 

4.1.  The  Multi-Channel  Detector  Experiment 

We  were  requested  to  conduct  experiments  to  determine  whether  the  multi-channel  detector 
can  be  enhanced  by  using  p-level  wavelet  bases  generated  by  p^-level  wavelet  functions  that  are 
/2-closest  to  Doppler  shifts  or  multi-bounced  versions  of  the  prototype  signal.  We  ran  the  multi¬ 
channel  detector  using  different  wavelet  bases  to  see  if  the  appropriate  basis  will  have  a  better 
chance  of  detecting  the  expected  return  signal,  compared  to  other  wavelet  basis.  A  thorough 
discussion  of  the  experimental  results  is  given  in  the  following  chapter.  This  subsection  will 
explain  the  experimental  procedures.  The  results  of  the  tests  are  reported  in  the  following  chapter. 

The  first  step  was  to  construct  the  detection  basis  for  the  different  Doppler  shifts  of  the 
signal,  CHIP2.  To  do  this,  we  ran  the  signal  through  the  TSS,  assuming  only  Doppler  of  0  and  +5 
kts.  The  return  signals  were  then  used  to  generate  the  /2-closest  3-level  wavelet  bases.  This  forms 
the  standard  and  the  Doppler-only  detector. 

The  next  step  was  to  construct  the  detection  basis  to  match  a  multi-path  signal  (i.e.,  a  signal 
affected  by  bottom  bounce).  This  signal  is  constructed  by  sending  the  prototype  signal  through  the 
TSS  with  default  setting  on  every  module  except  the  multi-path.  The  parameters  for  the  multi-path 
are: 


•  Sea  Depth  =  2000  ft, 

•  Signal/Receiver  Source  Depth  =  1000  ft, 

•  Target  Depth  =  1000  ft, 

•  Source-to-Target  Range  =  30  nm. 

These  are  also  our  standard  settings  for  the  multi-path  module. 

We  proceed  as  with  the  Doppler-only  case  and  generate  the  /  -closest  3-level  wavelet  basis 
to  the  multi-path  signal. 
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The  third  step  was  to  built  a  detection  basis  to  match  a  Doppler-shifted,  multi-path  signal. 
The  Doppler  shift  is  +5  kts,  and  the  multi-path  parameters  are  the  same  as  above.  Again,  we 
generate  the  /^-closest  3-level  wavelet  basis  for  this  signal. 

We  now  run  the  experiment,  sending  a  signal  through  the  TSS,  with  and  without  Doppler, 
and  with  and  without  multi-path.  The  Doppler  range  for  the  experimental  signals  are  0  kts.  ±5  kts, 
±10  kts.  The  multi-path  is  set  to  the  above  parameters.  Recall  that  the  return  signal  consists  of  9 N 
=  9216  data  points.  The  multi-channel  detector  assumes  Doppler  speeds  of  0  kts,  ±2  kts,  ±4  kts, 
±6  kts,  ±8  kts,  and  ±10  kts  -  making  up  the  eleven  channels  of  the  detector. 

As  mentioned  in  the  above  section  on  the  multi-channel  detector,  the  output  of  the  detector 
is  8  N  =  8192  data  point  per  channel.  It  is  impossible  for  the  printer  to  get  the  resolution  required 
for  eleven  channels  of  8192  data  points,  so  we  resample  the  data  points  with  a  peak-picking 
method. 

The  first  step  is  to  set  a  threshold  so  that  all  data  points  below  the  threshold  are  set  to  zero. 
Next,  we  partition  the  8192  data  points  into  groups  of  10  (in  the  order  they  were  originally 
presented).  In  each  of  the  groups  of  10,  we  pick  the  largest  value  (i.e.,  peak-picked)  and  use  it  as 
the  representative  for  that  group  of  10.  Thus,  we  get  819  peak-picked  data  points.  These  points 
are  then  plotted  on  a  grass  plot.  A  grass  plot  is  a  plot  of  a  sequence,  cut  below  by  the  threshold. 

We  also  conducted  tests  to  show  the  efficacy  of  /2-closest  wavelet  bases  for  compressed 
Neyman-Pearson  detection.  We  compared  the  /2-closest  wavelet  basis  to  the  real  Shannon  wavelet 
basis,  the  Fourier  basis,  the  Daubechies'  D4  wavelet  basis,  and  the  delta  (time)  basis.  The  tests 
were  to  determine  how  much  dB  loss  were  incurred  by  compressing  the  detecting  basis  for  the 
various  bases.  We  ran  these  tests  for  the  CHIP2  signal  sent  through  the  TSS  with  the  following 
settings: 


•  Only  Doppler,  at  0  kts,  ±2  kts,  ±4  kts 

•  Only  Multi-Path,  at  standard  parameters 

•  Doppler  and  Multi-path, 

•  Translated  signal 
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4.2.  Compression  Loss  versus  Compression  Ratio  Tests 


To  implement  the  compressed  Neyman-Pearson  (N-P)  detector,  we  need  to  find  bases  (in 
the  sense  of  vector  space  bases)  that  do  not  degrade  much  due  to  compression,  indeed,  the.  notion 
of  a  basis  is  equivalent  to  the  detector.  We  show,  in  the  following  chapter,  that  the  /"-closest  (to 
the  prototype  signal  CHIP2)  wavelet  basis  significantly  outperforms  all  other  bases.  However,  we 
do  not  have  mathematical  proof  that  this  basis  is  optimal  for  compressed  detection.  Additionally, 
we  do  not  have  proof  that  the  compressed  N-P  detector  is  near  optimal  (see  Chapter  VIII,  Section 
3). 


The  N-P  detector  is  optimal  (see  Chapter  IV),  but  it  is  slow.  Compressing  the  data 
increases  the  computation  speed,  the  drawback  being  some  signal  loss.  This  can  be  measured  as  a 
dB  loss,  called  the  degradation  due  to  compression.  For  different  bases,  the  amount  of  degradation 
loss  differs  for  an  equal  amount  of  compression. 

The  bases  under  comparison  are  the  3-level  /2-closest  wavelet  basis,  the  3-level  real 
Shannon  wavelet  basis,  the  3-level  Daubechies  D4  wavelet  basis,  the  Fourier  basis,  and  the  delta 
basis.  The  Fourier  basis  is  actually  the  absolute  value  of  the  Fourier  transform.  The  /2-closest 
wavelet  basis  is  generated  from  the  3r(*-level  wavelet  function,  which  is  /2-closest  to  the  prototype 
signal,  CHIP2. 

We  test  these  bases  for  degradation  under  various  scenarios;  we  consider  detecting  the 
prototype,  translates  of  the  prototype,  and  Doppler  shifts  of  the  prototype.  Not  all  bases  will  be 
tested  under  all  these  scenarios.  For  example,  there  is  no  need  to  test  the  Fourier  basis  under 
translations  since  the  absolute  value  of  the  Fourier  transform  is  translation  invariant. 

To  quantitatively  measure  how  good  a  particular  basis  is  for  compressed  Neyman-Pearson, 
we  compute  the  dB  loss  due  to  compression.  Because  the  Neyman-Pearson  detector  is  basis- 
independent,  all  bases  will  give  optimal  detection  if  used  without  compression.  What  we  want  to 
measure,  then,  is  the  entropy  of  the  transform. 

Let  B  be  a  particular  orthonormal  basis  (e.g.,  the  /2-closest  wavelet  basis,  real  Shannon 
wavelet  basis)  into  which  we  want  to  transform  the  signal.  If 

s  =  s[n|  =  (s(0),  s(l) . s (N  -  1)) 
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is  the  signal  in  the  time  basis,  and 


S  =  £[n]  =  (s(0),  s(  1 ), s(N  -  1)) 

is  the  signal  in  basis  %  we  rearrange  the  absolute  value  of  the  coefficients. 

{  |s(0)|,  |s(  1  )|, ....  |s(/V  -  1)| }, 

in  decreasing  order  so  that  choosing  the  M  largest  weights  is  equivalent  to  taking  the  first  M 
coefficients;  these  are  then  the  elements  of  A.  Assume  now  that  the  coefficients  are  reordered. 
Using  the  equations  from  Chapter  IV,  we  compute  the  entropy: 

(M- 1  U 

ZA=  X  ls(kf  f 

\  k  =  0  / 

and  then  to  determine  a  dB  loss,  we  take  the  logarithm: 

Da  =  20  log10  ZA, 


which  is  now  measured  in  dBs. 

In  our  tests,  \v:  up  the  runs  so  that  for  a  given  compression  ratio  (M  /  N),  the  code 

automatically  finds  the  optimal  subset.  A,  of  cardinality  M.  The  output  is  then  a  dB  loss  for  every 
value  of  M.  from  1  to  N .  The  graph  of  dB  loss  versus  compression  ratio  is  a  monotonically 
increasing  function.  Because  of  the  semi-logarithmic  scale  of  the  axes,  the  graphs  do  not  look 
monotone  increasing  for  small  compression  ratios,  although  they  are. 

The  results  show  quite  clearly  that  the  /"-closest  wavelet  basis  stands  above  the  other  bases 
for  compressed  N-P  detection  purposes.  Chapter  VII  shows  these  results. 
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CHAPTER  VII 


SUMMARY  OF  THE  EXPERIMENTAL  RESULTS 

We  present  and  discuss  the  results  of  the  experiments  and  tests  described  in  the  previous 
chapter.  Our  main  conclusion  for  the  multi-channel  detector  experiment  is  that  the  results  are 
inconclusive.  The  conclusions  for  the  compression  loss-versus-compression  ratio  tests  are  more 
conclusive;  they  do  point  out  that  /  -closest  wavelet  bases  perform  better  than  other  orthonormal 
bases  with  regard  to  the  compressed  Neyman-Pearson  detector.  The  test  results  are  illustrated  by 
graphs.  The  multi-channel  detector  experiment  results  are  also  graphically  illustrated.  Section  1 
presents  the  results  of  the  multi-channel  detector  experiment,  and  Section  2  presents  the  results  of 
the  compression  loss-versus-compression  ratio  tests. 

I.  The  Multi-Channel  Detector  Experiment 

The  results  of  the  multi-channel  detector  experiment  are  inconclusive.  We  wanted  to  know 
whether  the  multi-channel  detector  could  be  enhanced  by  using  the  Doppler  shift  of  the  signal 
(CHIP2)  as  the  detecting  prototype  signal  (see  Chapter  VI,  Section  4.1  for  details).  The  results 
show  no  appreciable  improvement  using  the  /2-closest  (to  the  Doppler  shift  of  the  signal)  wavelet 
basis,  compared  to  the  /  -closest  (to  the  prototype  signal)  wavelet  basis.  However,  the 
inconclusive  result  should  not  be  misconstrued  as  that  the  detector  is  defective.  The  results  do 
show  that  the  detector  does  work.  This  is  evident  in  the  noticeable  spikes  in  the  grass  plots.  It  is 
just  that  our  intent  was  to  see  whether  the  detector  can  be  enhanced  by  choosing  a  more  appropriate 
detection  basis.  Moreover,  the  result  does  not  disprove  the  efficacy  of  tailoring  the  detecting  basis 
to  the  expected  return  signal. 

The  explanation  for  the  inconclusive  result  is  straightforward:  the  signal,  CHIP2,  is 
relatively  Doppler-invariant.  As  mentioned  in  Chapter  V,  HFM  signals  are  invariant  under 
Doppler  shifts.  Since  CHIP2  is  constructed  by  concatenating  a  series  of  HFM  signals,  it,  too.  is 
Doppler  invariant.  Compare  the  time-series  plots  of  the  prototype  signal  and  a  +5  kt  Doppler  shift 
of  the  signal  (see  Figure  VI-7).  The  two  plots  are  almost  identical.  The  time-series  of  the  shallow- 
channel  multi-path  effects  (with  and  without  +5  kt  Doppler)  are  also  given  in  Figure  VI-7;  note  the 
differences  between  the  signals  with  and  without  multi-path. 
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We  also  plot  the  T-F  distributions  of  these  four  return  signals,  using  the  windowed  Fourier 
Spectrogram,  as  described  in  Section  3  of  Chapter  V  (see  Figures  VII- 1,  VII-2,  VII-3,  and  VIM). 
The  distributions  verify  again  the  similarities  of  the  Doppler  shifts  to  the  unshifted  versions. 

The  results  of  the  multi-channel  detector  are  read  from  the  grass  plots.  Figures  VII-5,  VII- 
6,  VII-7,  and  VII-8  show  some  grass  plots  from  the  experiment.  These  four  plots  all  have  a 
compression  ratio  equal  to  0.0 1 .  For  clarity,  the  grass  plots  do  not  include  the  first  window,  which 
is  comparing  noise,  nor  the  fifth  through  eighth  windows,  which  include  the  low  end  of  the  signal 
in  addition  to  noise.  We  have  test  results  using  various  possible  combinations  of  detector  signal 
(CHIP2,  CHIP2  with  Doppler,  CHIP2  with  multi-path,  CHIP2  with  Doppler  and  multi-path), 
return  signal  (similar  to  detector  signal  choices),  and  number  of  weights,  and  we  also  compared 
CHIP2  results  with  a  pulse  result.  For  brevity,  we  decided  not  to  include  all  these  results. 

In  the  grass  plots  (those  from  the  figures  shown  in  this  report  and  others),  there  is  a 
noticeable  peak  corresponding  to  the  correct  matching  of  the  return  signal  with  the  detecting  signal. 
Some  of  the  plots  (mainly  those  in  which  the  detecting  basis  is  generated  from  the  multi-path 
version  of  the  signal  or  the  return  signal  has  been  sent  through  the  multi-path  module)  show  that 
the  main  peak  is  not  as  definite.  This  suggests  that  the  basis  choice  (hence,  the  signal  choice)  is  not 
a  good  choice.  It  is  in  the  similarity  between  the  various  channels  that  demonstrates  the 
inconclusiveness  of  the  whole  experiment. 

Let  us  study  the  plots  in  Figures  VII-5  through  VII-7.  Figure  VII-5  shows  the  result  of 
attempting  to  detect  a  non-Doppler-shifted  single-path  return  CHIP2  signal  with  a  detector  basis 
built  from  the  CHIP2  signal.  The  peak  occurs  at  the  beginning  of  the  third  window,  which  is  the 
first  "signal"  window.  Note  the  similarities  between  the  different  channels.  This  is  due  to  the  fact 
that  CHIP2  is  Doppler-invariant;  hence,  the  detector  bases  for  the  channels  will  be  very  similar  to 
each  other.  Indeed,  in  all  of  the  four  figures  (and  many  others  not  illustrated  in  this  report) 
variations  between  channels  are  minimal,  due  to  the  Doppler  invariance  of  CHIP2. 

The  major  differences  occur  when  comparing  different  detector  bases,  most  noticeably 
when  we  use  a  CHIP2  with  multi-path  signal  to  generate  the  detector  basis.  The  grass  plot  of 
Figure  VII-6  (which  uses  CHIP2  with  multi-path  to  generate  the  detector  basis)  starkly  contrasts 
with  that  of  Figure  VII-5.  The  lack  of  a  noticeable  peak  in  Figure  VII-6  shows  that  such  a  basis  is 
not  effective. 
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Unfortunately,  the  detector  basis  generated  by  CHIP2  with  multi-path  is  equally  ineffective 
in  attempting  to  detect  a  CHIP2  with  multi-path  return  signal.  This  is  seen  in  Figure  Vll-7.  Note 
again  the  lack  of  a  distinctive  peak.  Our  results  indicate  that  the  /2-closest  wavelet  basis  built  out  of 
CHIP2  with  multi-path  is  not  a  good  detector  basis. 

Figure  VII-8  shows  a  prominent  peak,  despite  attempting  to  detect  a  non-Doppler-shifted 
signal.  Again,  this  illustrates  the  Doppler  invariance  of  CHIP2. 

In  general,  our  results  show  that  detecting  multi-path  using  a  detector  basis  built  from 
either  CHIP2.  CHIP2  with  Doppler,  or  CHIP2  with  multi-path  was  ineffective.  Not  one  of  the 
plots  show  any  prominent  peaks,  and  even  increasing  the  number  of  weights  did  not  improve 
detection  discrimination.  Compare  Figures  Vll-7,  VII-9,  and  VII- 10.  This  lack  of  distinction  also 
occurs  when  CHIP2  with  multi-path  is  used  to  build  the  detector  basis.  Figures  VII- 1 1  and  VII- 
12  illustrate  these  grass  plots. 

We  attempted  using  the  multi-channel  detector  on  other  signals.  One  other  signal  choice  is 
the  PULSE  100  (short)  signal.  The  results  were  also  not  very  encouraging  (see  Figure  VII- 13. 
The  Doppler  Ambiguity  Displays  of  pulse  signals  (see  Figures  V-17  and  V-18)  show  that  the 
effects  of  Doppler  are  appreciable  only  for  shifts  corresponding  to  speeds  greater  than  +10  kts. 
Since  our  multi-channel  detector  output  for  the  PULSE100  (short)  is  limited  between  -10  kts  to 
+  10  kts,  this  may  be  the  reason  for  the  lack  of  visual  differences. 

Our  conclusion  is  that  an  appropriate  choice  for  the  detector  basis  is  not  necessarily  the 
/2-closest  wavelet  basis  to  the  expected  return  signal.  It  is  also  possible  that  certain  return  signals 
will  be  difficult  to  detect,  no  matter  what  the  detector  basis  is.  We  did  not  try  using  more  levels  in 
our  detector  basis,  which  might  have  enhanced  detection.  (This  may  be  one  avenue  of  short-term 
future  research.)  The  multi-channel  detector  experiments  were  all  conducted  using  3-level  wavelet 
bases.  We  know  higher-level  wavelet  bases  yield  better  compressed  detection  capabilities.  At  this 
stage,  we  do  not  know  whether  the  cause  for  the  poor  performance  is  in  the  choice  for  the  detection 
basis  or  in  something  more  fundamental.  We  feel  the  accuracy  of  the  simulation  system  (TSS)  is 
sufficient  to  warrant  dismissing  it  from  consideration  as  a  possible  cause  for  the  poor  performance 
of  the  multi-channel  detector. 
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Fourier  Spectragram  of  CHIP2 
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Fourier  Spectragram  of  +5  kts  Doppler  Shift  of  CHIP2 
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Fourier  Spectragram  of  Multi-Path  of  CHIP2 
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Fourier  Spectragram  of  +5  kts  Doppler  Shift  and  Multi-Path  of  CHIP2 
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Compressed  Detector  +10  to -10  Doppler  10  Weights  with  Peak  Picking 
10*log10(abs(SPL))  Range:  200dB  to  220dB  9. 10  seconds  at  1000  samples  /  second 

Actual:  0  Doppler  Single  Path 


Compressed  Detector  +10  to -10  Doppler  10  Weights  with  Peak  Picking 

10*log10(abs(SPL))  Range  200dB  to  220dB  9.18  seconds  at  1000  samples  /  second 
Actual:  Chip2/Multi-Path  0  Doppler  Single  Path 
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Actual:  Chip2/Doppler  0  Doppler  Single  Path 
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Actual:  Chip2/Multi-Palli  +  5  Doppler  Single  Path 
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2.  Compression  Loss-versus-Compression  Ratio  Tests 

The  results  of  these  compression  loss-versus-compression  ratio  tests  are  much  more 
encouraging  and  substantial  than  the  multi-channel  detector  experiment.  In  these  series  of  tests,  we 
have  definitely  determined  that  the  /"-closest  wavelet  basis  outperforms  all  other  standard  bases  for 
compressed  Neyman- Pearson  detection. 

Recall  from  the  previous  chapter  that  the  goal  of  choosing  an  appropriate  basis  is  borne  out 
by  little  compression  loss  for  large  compression  ratios.  Figure  VII- 14  shows  a  hypothetical 
compression  loss-versus-compression  ratio  plot,  with  directions  showing  what  makes  for  a 
"good"  or  "bad"  detection  basis. 


Figure  VII- 14 


Compression  Loss  (dB) 


In  using  the  basic  compressed  Neyman-Pearson  detector,  we  do  not  zero-pad  the  signal; 
hence,  we  create  the  /"-closest  3-level  wavelet  basis  based  on  the  signal.  See  Figure  VI-4  for  the 
time-series  plot  of  the  l2- closest  3r^-level  wavelet  function.  Figure  VI-5  shows  the  difference 
between  the  CHIP2  signal  and  the  /"-closest  3rc^-level  wavelet  function.  Figure  VI-6  shows  the 
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time-series  plots  of  all  the  wavelet  functions:  first,  second,  third,  and  the  extra  third-level  wavelet 
functions. 


Figure  VII- 15  shows  the  wavelet  transform  of  CHIP2  (sometimes  abbreviated  as  CH2), 
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offset  by  0,  1,2,  and  3.  The  wavelet  basis  is  the  / -closest  3-level  wavelet  basis.  One  can  see  that 
the  signal  is  highly  concentrated  in  the  zero  offset  at  n  =  832.  We  also  sent  this  signal  through  our 
transmission  simulation  system,  simulating  a  +5  kt  Doppler,  a  shallow-channel  bottom-bounce 
effect,  and  the  combination  of  both:  a  +5  kt  Doppler  shift  and  the  effect  of  a  shallow  channel. 
Figure  VI-7  shows  the  time-series  plots  of  CHIP2  in  these  situations,  compared  to  the  prototype. 
Figure  VII- 16  shows  the  wavelet  transform  of  each  of  the  above  affected  (Doppler,  shallow 
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channel)  signal,  with  respect  to  the  /‘-closest  3-level  wavelet  basis.  As  expected,  the  prototype 
signal  shows  the  best  concentration  of  the  weights  in  a  small  interval,  and  the  shallow  channel  and 
+5  kt  Doppler  case  has  tl  e  worst  concentration  of  weights. 

For  comparison  to  the  other  bases.  Figure  VII- 17  shows  the  transforms  of  CHIP2  in  the 
delta  basis  (time  basis,  hence  the  identity  transform),  the  absolute  part  of  the  Fourier  transform, 
and  the  real  Shannon  wavelet  transform,  and  the  /2-closest  wavelet  basis.  As  expected,  the 
/2-closest  wavelet  basis  shows  the  best  concentration  of  the  weights  in  a  small  interval. 

Figures  VII- 18  through  VII-33  show  the  compression  loss  curves  for  various  bases  for  the 
purpose  of  comparison.  Figure  VII- 18  compares  five  bases  in  the  case  when  there  is  no  Doppler- 
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shifting  and  the  signal  is  not  offset  in  the  window;  note  that  the  /  -closest  wavelet  basis  is  the  clear 
winner.  Figures  VII- 19  through  VII-23  show  the  result  of  offsetting  the  signal  in  the  window  on 
the  various  individual  compressed  detectors.  Figure  VII-24  compares  the  /2-closest  wavelet 
basis,  the  D4  wavelet  basis,  and  the  real  Shannon  wavelet  basis  when  the  signal  is  offset  by  2. 
Figure  VII-25  compares  the  /2-closest  wavelet  basis,  the  real  Shannon  wavelet  basis,  and  the 
Fourier  basis  when  all  offsets  from  0-4  have  been  combined  to  form  worst-case  curves.  (The 
recursive  nature  of  the  wavelet  basis  implies  only  the  0  through  4  translates  be  studied;  any  greater 
translates  will  yield  redundant  results.)  In  both  Figures  VII-25  and  VII-25,  the  /2-closest  wavelet 
basis  performs  the  best.  Figures  VII-26  through  VII-3I  show  the  result  of  Doppler-shifting  on 
six  different  bases.  Figure  VII-32  compares  the  /2-closest  wavelet  basis,  the  D4  wavelet  basis,  the 
delta  basis,  the  Fourier  basis,  and  the  real  Shannon  basis  when  there  is  Doppler-shifting 
corresponding  to  -4  knots.  Figure  VII-33  compares  the  l2- closest  wavelet  basis,  the  delta  basis, 
the  Fourier  basis,  and  the  real  Shannon  basis  when  Doppler  shifts  from  -4  knots  to  +4  knots  (in  2- 
knot  increments)  have  been  combined  to  form  worst-case  curves.  In  these  last  two  figures,  the 
/2-closest  wavelet  basis  is  once  again  the  winner. 
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Finally,  the  most  important  test  of  computation  speed  is  the  amount  of  real  time  it  takes  to 
perform  the  computations  to  achieve  a  certain  level  of  degradation.  We  conducted  timing 
comparison  studies  and,  again,  showed  the  /“-closest  wavelet  basis  achieves  a  low  amount  of 
degradation  in  a  shorter  computing  time  than  other  bases.  Figure  VII-34  shows  the  compression 
loss  as  a  function  of  relative  time. 
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CHAPTER  VIII 


CONCLUSIONS  FROM  THE 
PRESENT  RESEARCH 


We  summarize  the  main  results  achieved  in  our  research  program  on  LFA/LPI  signals. 
This  chapter  is  divided  into  two  main  sections,  the  first  concerns  the  experimental  results  and  the 
second  concerns  the  theoretical  results.  We  will  discuss  improvements  and  future  wo.k  in  the 
following  chapter. 

1.  Experimental  Results 

Our  efforts  have  generated  the  following  experimental  results: 

•  A  compressed  detection  method  and  a  detector  that  suffers  the  least  degradation  due  to 
compression,  in  head-to-head  comparison  to  other  detector  choices. 

•  A  multi-channel  detector  for  simultaneously  detecting  possible  Doppler  and  time  shifts 
of  the  return  signal,  relative  to  the  prototype  signal. 

•  A  transmission  simulation  system  (TSS)  that  simulates  all  phases  and  major 
environmental  effects  encountered  by  underwater  sonar  signals  during  their 
propagation. 

Compressed  Detector.  Our  compressed  detector  is  based  on  the  compressed  Neyman- 
Pearson  detector  (see  following  section  on  theoretical  results).  We  tested  several  bases  to 
determine  which  performed  better.  The  results  (see  Chapter  VI  for  details  on  the  tests  and  Chapter 
VII  for  the  results)  show  the  l2 -closest  wavelet  basis  significantly  outperforms  all  other  bases, 
including  other  standard  wavelet  bases,  like  the  Daubechies  D4  basis  and  the  real  Shannon  basis. 
The  /’-closest  wavelet  basis  detector  is  essentially  a  matched  filtering  system.  It  does  have 
drawbacks,  which  we  acknowledged  and  dealt  with  through  the  use  of  the  binary  decomposition 
method  for  detection  purposes. 

We  emphasize  that  our  compressed  detector  is  indeed  compressed:  of  the  1024  data  joints 
describing  a  finite  discrete  signal  picked  up  by  the  detecting  window,  we  use  on  the  order  of  to 
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16  data  points  for  detection  purposes  and  have  losses  that  are  less  than  2  -  3dB  overall!  Figures 
VII-(  18)  though  VII-(33)  show  that  the  detector  using  the  /2-closest  wavelet  basis  suffers  very 
little  degradation  compared  to  detectors  using  the  Fourier,  the  delta,  or  other  wavelet  bases.  In 
order  to  perform  compression,  the  signal  must  be  transformed  into  the  selected  wavelet  basis.  We 
have  shown  that  this  can  be  done  rapidly  with  little  impact  on  the  speed  gains  due  to  compression. 

Simulation.  The  transmission  simulation  system  allowed  us  to  simulate  Doppler  and 
time-delay  effects  on  prototype  signals.  We  also  wanted  the  prototype  signals  to  be  relevant  to 
actual  sonar  signals  and  so,  we  developed  a  program  to  generate  these  signals.  The  combination 
of  having  a  Doppler  simulator  and  real-life  signals  allowed  us  to  test  the  effectiveness  of  the 
detector.  Our  research  in  this  direction  led  us  to  develop  the  multi-channel  detector,  which  is 
capable  of  simultaneously  detecting  possible  Doppler  and  time  shifts  of  the  signal. 

The  transmission  simulation  system  involves  over  10,000  source  lines  of  code.  Our 
updated  version  includes  the  Pecris  Model  to  deal  with  shallow  water  channels,  in  which  bottom 
bounces  have  an  appreciable  effect  on  the  signal.  The  simulation  is  also  equipped  to  compute 
Doppler  ambiguity  displays  as  well  as  wavelet  and  Fourier  waterfall  displays.  The  transmission 
simulation  system  is  described  in  Chapter  VI.  It  figures  prominently  in  testing  the  multi-channel 
detector  and  other  detection  methods  (see  Chapter  VII  for  the  results  of  these  tests). 

Multi-channel  Detector.  The  multi-channel  detector  performed  as  intended,  emphasis  on 
the  recursive  nature  of  wavelets  led  to  fast  computation  and  minimal  storage  requirements. 
Unfortunately,  experiments  conducted  with  the  multi-channel  detector  did  not  illustrate  its  full 
capabilities,  due  to  the  poor  choice  of  prototype  signals.  We  address  future  experiments  on  the 
multi-channel  detector  in  the  next  chapter. 

2.  Theoretical  Results 

In  the  theoretical  arena,  we  have  researched  into  the  following: 

•  Near-optimal  detectors,  based  on  the  optimal  Neyman-Pearson  detector. 

♦  A  binary  decomposition  method  for  signal  detection  to  address  the  signal  translation 
issue. 
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•  New  approaches  to  using  discrete  wavelets  for  detection  methods  as  well  as  for 
perfecting  time-frequency  distributions. 

•  A  better  understanding  of  signals  from  using  time-frequency  analysis. 

Detection.  During  our  research  efforts,  we  originally  devised  a  detection  method  which 
was  eventually  shown  to  be  suboptimal.  Our  research  showed  the  Neyman-Pearson  detection 
method  is  optimal,  and  thus,  should  be  implemented.  However,  the  Neyman-Pearson  detector 
requires  enormous  computing  efforts  not  available  for  real-time  detection  and  immediate 
decisions.  We  then  developed  the  compressed  version  of  the  detector  so  that  this  compressed 
version  became  much  faster. 

In  conjunction  with  the  detector  is  the  detector  basis.  The  choice  of  basis  may  allow  for 
additional  compression  without  significant  degradation.  Through  wavelets,  we  discovered  that  the 
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/  -closest  wavelet  basis  yields  very  little  loss  even  at  high  compression. 

Binary  Decomposition.  Other  uses  of  wavelets  include  the  binary  decomposition 
method,  which  has  the  additional  capability  of  detecting  translates  of  the  signal  better  than  just  the 
/2-closest  wavelet  basis.  We  also  studied  the  possibility  of  creating  noise-optimal  wavelets.  Such 
wavelet  bases  would  discriminate  the  signal  from  the  noise  by  "matching"  the  wavelets  to  the 
noise.  Then,  the  wavelet  coefficients  corresponding  to  the  noise  would  be  concentrated  over  a 
small  support,  disjoint  (hopefully)  from  the  support  for  the  signal  coefficients.  Our  results  here 
showed  that  such  bases  were  highly  dependent  on  the  noise  and  often  affected  the  signal  negatively 
and  thus,  such  wavelets  were  eventually  discarded. 

Our  work  in  developing  signals  created  by  wavelets  showed  that  this  procedure  was  indeed 
possible,  but  the  process  does  not  match  the  capabilities  of  current  sonar  signal  generation 
methods.  Recently,  we  have  developed  a  wavelet  spectrogram  from  our  work  in  time-frequency 
analysis.  This  new  tool  may  revive  our  signals  generation  goal  with  a  new  procedure  that  will 
definitely  be  fast  and  can  be  tailored  to  have  specific  characteristics.  We  will  elaborate  upon  this  in 
the  following  chapter. 

Time-Frequency  Analysis.  Time-frequency  analysis,  performed  by  computing  a  time- 
frequency  distribution,  is  a  viable  tool  for  investigating  the  properties  of  a  signal,  especially  highly 
non-stationary  signals,  such  as  transient  sonar  signals.  Our  research  into  this  area  came  only 
recently.  We  plan  for  further  work  in  this  area,  as  described  in  the  next  chapter. 
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CHAPTER  IX 


RECOMMENDATIONS  FOR  FUTURE  RESEARCH 


In  this  chapter,  we  discuss  areas  in  which  further  work  will  benefit  DARPA's  LFA 
program. 

1.  Experimental 

This  area  addresses  testing  of  the  compressed  wavelet  detector  and  the  multi-channel 
detector.  To  date  all  of  the  testing  of  the  detectors  has  been  with  simple  simulation  models  such  as 
the  Pecris  model.  Further  testing  is  necessary  to  demonstrate  the  effectiveness  of  the  detectors  in  a 
real-world  scenario.  The  following  is  recommended. 

•  Incorporate  the  Safari  Model  into  the  LFA  Simulation. 

•  Perform  testing  of  the  compressed  detector  and  the  multi-channel  detector  using  the 
Safari  model  to  simulate  the  sound  channel 

•  Perform  testing  of  the  compressed  detector  and  the  multi-channel  detector  using  a 
sound  channel  derived  from  a  real  world  experiments  using  an  impulse  charge  at- 
sea. 

•  Perform  test  of  the  compressed  detector  and  the  multi-channel  detector  using  the 
derived  real  world  channel  and  signals  actually  generated  and  received  at-sea. 

•  Determine  the  sensitivity  of  the  multi-channel  detector  to  Doppler  shifts  using  a 
signal  with  greater  Doppler  dependence. 

For  example,  our  use  of  wavelets  over  the  span  of  this  research  phase  has  evolved  from 
signal  generation  to  signal  processing.  We  now  feel  this  is  the  proper  role  for  wavelets  in  regard  to 
LFA/LPI  signals.  We  also  have  reason  to  believe  that  wavelets  will  provide  a  new  means  to 
signal  generation,  one  that  we  did  not  consider  during  our  current  research  phase. 
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2.  New  Development 


Future  theoretical  research  topics  include: 

•  Signal  generation  using  the  wavelet  spectrogram.  This  effort  is  long-term. 

•  Time-frequency  distribution  for  non-stochastic  noise  analysis.  This  work  is  definitely 
short-term. 

•  Additional  utilizations  of  the  /2-closest  wavelet  basis.  This  study  is  potentially  a  long¬ 
term  research  problem. 

•  Other  detection  methods.  This  research  topic  is  most  likely  a  long-term  research 
proposal. 

•  Technical  merits  of  the  binary  decomposition  method.  This  effort  is  short-term. 

We  were  asked  to  study  the  possibility  of  creating  signals  having  prescribed  Doppler 
Ambiguity  displays  (see  Chapter  V,  Section  2  for  a  brief  description  of  Doppler  Ambiguity  plots). 
We  were  able  to  construct  signals  having  certain  properties,  as  suggested  by  their  Doppler 
Ambiguity  plots,  but  we  have  not  been  able  to  take  an  ambiguity  plot  and  create  a  signal  from  it. 
We  have  made  some  progress  in  this  area  and  see  this  as  a  possible  short-term  effort. 

Although  our  use  of  wavelets  has  moved  from  signal  generation  to  signal  processing,  we 
do  see  a  possibility  of  the  application  of  wavelets  in  signal  generation.  This  application  relates  to 
the  wavelet  spectrogram.  One  theoretical  advantage  of  the  wavelet  spectrogram  is  that  it  is  an 
invertible  process.  While  all  the  other  T-F  distributions  are  non-invertible  processes  (i.e.,  many 
signals  may  yield  the  same  T-F  distribution  plots,  due  to  taking  the  absolute  value  and  due  to  the 
localizing  nature  of  the  distribution  functions),  the  wavelet  spectrogram  is  an  invertible  process; 
for  each  signal,  there  is  an  unique  wavelet  spectrogram  plot,  and  conversely.  This  means  there  is  a 
quick  way  to  construct  a  signal  with  a  prescribed  wavelet  spectrogram.  More  importantly,  one 
must  understand  that  this  process  cannot  be  accomplished  with  any  of  the  other  T-F  distributions. 
In  this  aspect,  signal  generation  is  a  long-term  effort. 

This  use  of  wavelets  will  enable  us  to  model  noise  more  effectively  than  our  current 
moving-average  methods.  Certain  noises,  like  biologies  and  reverberation,  are  not  completely 
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stochastic;  there  are  certain  key  frequencies  that  characterize  the  noise.  These  frequencies  can  be 
determined  through  T-F  distribution  plots,  and  time-series  noise  can  be  recreated  by  an  inverse 
wavelet  spectrogram. 

For  short-term  efforts,  we  would  like  to  apply  the  detector  on  the  T-F  distribution  data. 
Moreover,  we  have  worked  on  a  classification  method  to  classify  families  of  signals.  We  would 
like  to  apply  this  classification  method  to  T-F  distributions.  In  this  way,  we  can  classify  signal 
types  by  the  signal's  spectrogram. 

Finally,  we  want  to  determine  some  of  the  technical  merits  of  the  \j/-binary  decomposition 
algorithm  for  detection  purposes.  The  y-binary  decomposition  algorithm  was  developed  late  in 
the  research  phase  and  has  not  been  fully  tested.  Testing  it  will  be  a  short-term  effort. 

We  would  like  to  consider  other  detection  methods.  Although  the  Neyman-Pearson 
detector  is  optimal,  the  compressed  Neyman-Pearson  is  only  near-optimal.  Perhaps  there  is  an 
alternative  detection  method  that  yields  better  results  for  equal  or  comparable  computing  effort. 

3.  Conclusion 

We  believe  that  our  research  to  date  has  demonstrated  the  utility  of  wavelets  in  the  fields  of 
signal  generation,  signal  processing,  detection,  and  classification.  We  have  developed,  in 
particular,  detection  algorithms  that  have  greater  accuracy  and  greater  speed  than  the  standard 
techniques  in  use.  It  has  become  evident  from  our  research  on  this  project  that  wavelet  theory  is  an 
immensely  promising  and  fruitful  field  of  research  with  many  potential  applications;  continued 
research  is  truly  imperative.  It  is  our  hope  to  be  a  part  of  that  continued  effort  and  to  follow  up  on 
the  potential  topics  for  further  study  that  arose  during  the  research  performed  for  this  project. 
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